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ABSTRACT 


The  theory  of  oceanic  convection  and  entrainment  has  been  developed  mainly  in  horizontally 
homogeneous  regimes,  yet  large-scale  spatial  variability  is  known  to  control  the  sites  and  intensity  of 
deep  convection.  Wintertime  Greenland  Sea  conditions  were  selected  to  simulate  convection  and 
quantify  the  interplay  between  local  forcing  and  large-scale  gradients.  Here  circulation  and 
preconditioning  produce  horizontal  gradients  in  the  stratification;  some  of  the  resulting  stratification 
is  conducive  to  the  formation  of  thermobaric  convective  instabilities. 

A  large  eddy  simulation  (LES)  model  modified  to  include  large-scale  horizontal  density 
gradients  was  used  to  study  the  effects  of  the  gradients  on  turbulence.  Horizontal  turbulent  kinetic 
energy  (TKE)  and  scalar  variances  increased  compared  to  simulations  with  no  large-scale  gradient. 
The  additional  horizontal  TKE  is  created  at  scales  larger  than  the  convective  plume  scale.  A  mean 
horizontal  circulation  develops  in  response  to  the  large-scale  overturning.  The  balance  between 
convection  and  overturning  increases  stratification  in  the  lower  region  of  the  mixed  layer,  and  plumes 
may  undergo  slantwise  convection. 

One-dimensional  bulk  TKE  model  results  were  compared  to  a  large  eddy  simulation  of 
wintertime  Greenland  Sea  convection.  One-dimensional  and  LES  results  were  similar  in  the 
distribution  of  TKE  components  and  in  the  ratio  AR  of  entrainment  buoyancy  flux  to  surface  buoyancy 
flux  for  the  winter  period  modeled.  The  value  of  AR  was  large  because  of  strong  shear  production, 
0.42  for  the  one-dimensional  model,  and  0.36  for  the  LES.  Detraining  thermobaric  plumes  were 
simulated  by  LES  under  various  conditions  of  rotation  and  stratification.  A  critical  depth  h„  and 
critical  velocity  wcr  hypothesized  by  Garwood  et  al.  (1994)  were  shown  to  be  predictors  for  onset  of 
detrainment.  The  skewness  of  vertical  velocity  in  a  horizontal  slice  just  below  the  mixed  layer  is 
shown  to  be  an  indicator  for  detrainment  events.  The  fraction  of  mixed-layer  water  present  at  depth 
quantifies  plume  transport  below  the  layer. 
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I.  INTRODUCTION 


A.  STATEMENT  OF  PROBLEM 

All  interaction  of  the  atmosphere  with  the  ocean  must  begin  where  the  two  fluids  are 
in  contact  -  that  is,  in  the  adjacent,  fully  turbulent  layers  of  the  ocean  and  atmosphere.  The 
oceanic  mixed  layer  is  more  buoyant  than  the  water  below,  so  that  changes  in  the  atmosphere 
are  to  some  extent  isolated  from  the  deep  ocean  by  the  more  buoyant  layer.  Exceptions  to 
this  general  structure  are  in  the  polar  ocean  regions,  where  mixing  may  occur  to  great  depths, 
or  even  to  the  bottom.  This  coupling  of  the  ocean  to  the  atmosphere  through  deep  mixed 
layers  in  polar  regions  drives  the  global  thermohaline  conveyor  belt  and  thus  influences  the 
earth’s  climate  system. 

One-dimensional  models  are  widely  used  to  predict  mixed-layer  properties,  as  well 
as  the  depth  of  mixing.  The  depth  of  mixing  changes  due  to  entrainment,  that  is,  by  mixing 
into  it  the  nonturbulent  fluid  from  below.  Generally,  a  model’s  constants  are  tuned  to  give 
correct  entrainment  for  some  set  of  typical  mid-latitude  situations  and  will  not  be  as 
applicable  to  polar  and  subpolar  scenarios.  To  extend  the  range  of  conditions  for  which  a 
model  is  useful,  there  needs  to  be  comparison  to  data  from  a  wide  range  of  conditions.  Given 
the  scarcity  of  detailed  data  sets  of  turbulence,  especially  in  the  near-polar  regions,  the 
realistic  physics  of  large  eddy  simulation  (LES)  models  provides  important  data  for 
comparison.  Energy  distribution  among  the  components  and  the  relative  importance  of 
energy  sources  need  to  be  compared  between  one-dimensional  and  LES  models. 
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Nonlinearities  in  the  equation  of  state  for  seawater  add  complexity  to  mixed-layer 
dynamics  at  high  latitudes.  Near  freezing  temperature,  the  thermal  expansion  coefficient  is 
especially  sensitive  to  pressure;  this  effect,  termed  thermobaricity,  can  produce  instabilities 
when  a  cold  fresh  mixed  layer  overlies  warmer,  more  saline  water,  resulting  in  detrainment 
of  mixed-layer  water  into  the  ocean  interior.  Detraining  plumes  have  been  simulated  by 
Garwood  et  al.  (1994)  in  a  regime  with  neutral  stratification  below  the  mixed  layer.  It  has 
not  yet  been  shown  how  nonzero  stratification  will  affect  plume  formation,  plume 
entrainment,  and  the  interior  fluid. 

One-dimensional  models  and,  typically,  LES  simulations  of  ocean  turbulence  use  an 
assumption  of  large-scale  horizontal  homogeneity.  Yet  convection,  especially  the  selection 
of  sites  for  deep  convection,  is  known  to  be  affected  by  the  large-scale  ocean  conditions. 
There  is  a  need  for  basic  knowledge  of  how  horizontal  inhomogeneities  may  interact  with 
turbulence.  Energy  component  balance  and  energy  sources  for  convection  may  be  altered 
by  a  horizontal  gradient  or  by  vertical  current  shear,  and  stability  may  differ  from  the  stability 
found  in  horizontally  homogeneous  domains. 

B.  THE  GREENLAND  SEA  REGION 


Most  of  the  volume  of  the  world  ocean  is  insulated  from  the  atmosphere  by  the 
buoyant  upper  layer  of  the  ocean,  with  only  a  few  sites  subject  to  deep  mixing.  The  known 
sites  of  deep  convection  are  located  in  the  Weddell  Sea,  the  Meditteranean  Sea,  the  Labrador 
Sea,  and  the  Greenland  Sea.  As  early  as  1906,  Nansen  observed  extremely  small  differences 
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in  properties  of  Greenland  Sea  surface  and  deep  waters,  and  proposed  that  deep  water 
formation  occurred  there.  A  schematic  of  the  Greenland  Sea  is  shown  in  Figure  1.  In  the 
upper  layers,  warm,  salty  water  from  the  Norwegian  Atlantic  Current  is  advected  in  by  the 
West  Spitzbergen  Current.  Some  of  this  flow  is  recirculated  in  the  Return  Atlantic  Current, 
where  it  joins  the  cold,  fresh  East  Greenland  Current  flowing  in  from  the  Arctic  Ocean.  The 
eastern  edge  of  the  East  Greenland  Current  is  called  the  East  Greenland  Polar  Front.  The 
circulation  is  completed  with  the  cold  and  fresh  Jan  Mayen  Current.  The  cyclonic  circulation 
results  in  a  doming  of  isopycnals  in  the  central  Greenland  Sea  gyre.  This  brings  the  weakly 
stratified  interior  waters  close  to  the  surface  (Rudels  and  Quadfasel,  1991;  Schott  et  al., 
1993;  Worcester  et  al.,  1993)  and  is  the  first  step  in  preconditioning  for  deep  convection. 

In  early  winter,  as  the  marginal  ice  zone  propagates  rapidly  eastward,  brine  release 
and  entrainment  act  to  increase  the  salinity  of  the  surface  waters.  In  combination  with  strong 
surface  cooling,  the  buoyancy  of  the  mixed  layer  is  thus  eroded.  Mixed-layer  deepening 
during  this  period  may  extend  far  enough  into  the  warm  Atlantic  Intermediate  Water  so  that, 
when  subsequent  cooling  and  stirring  occur,  warming  by  entrainment  dominates  over  surface 
cooling  (Visbeck  et  al.,  1995).  Any  ice  present  is  rapidly  melted,  and  further  ice  formation 
is  prevented;  an  ice-free  area  known  as  the  Nordbukta  forms.  The  Nordbukta  is  transitory, 
and  is  bounded  on  the  west  by  the  East  Greenland  Polar  Front,  and  on  the  south  and  east  by 
a  cold  icy  tongue  called  the  Is  Odden.  The  Nordbukta,  situated  over  the  Greenland  Sea  gyre, 
is  further  preconditioned  for  deep  convection.  Within  the  preconditioned  region,  convection 
is  believed  to  be  localized  by  local  forcing,  migrating  cyclonic  or  anticyclonic  eddies,  or 
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other  unknown  mechanisms  (Hakkinen,  1987;  Hakkinen  et  al.,  1992;  Hakkinen,  1995; 
Kill  worth,  1979;  Killworth,  1983;  Lherminier,  1999;  Morawitz  et  al.,  1996). 

C.  OVERVIEW 


Chapter  I  introduces  the  topic  of  study  and  motivation  and  reviews  the  literature. 
Chapter  II  lays  out  the  fundamental  equations  and  modeling  techniques  used  in  this  study. 
Chapter  III  presents  the  case  study  of  a  wintertime  Greenland  Sea  scenario;  one-dimensional 
model  results  are  compared  to  an  LES  simulation  of  the  same  period.  Chapter  IV  presents 
LES  numerical  experiments  in  which  thermobaric  plumes  were  simulated  under  various 
conditions.  The  detrainment  of  plumes  from  the  mixed  layer  and  their  subsequent  mixing 
with  the  interior  are  affected  by  the  differing  interior  stratifications  used  in  several  model 
scenarios.  Chapter  V  shows  results  for  LES  simulations  of  active  convection  and  turbulent 
spin-down  in  the  presence  of  large-scale  horizontal  density  gradients.  Behavior  of  turbulent 
kinetic  energy  (TKE)  components  and  scalar  variances  differs  from  convection  in  a 
horizontally  homogeneous  regime.  Chapter  VI  summarizes  important  conclusions  and  makes 
recommendations  for  further  study. 
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D.  LITERATURE  REVIEW 


1.  Mixed  Layer  Theory  and  Modeling 

The  importance  of  ocean  mixed-layer  physics  is  expressed  by  Garwood  (1976): 

Interest  in  the  ocean  mixed  layer  stems  from  both  theoretical  and 
practical  considerations.  Thermal  energy  and  mechanical  energy  received 
from  the  atmosphere  not  only  control  the  local  dynamics,  but  the  layer  itself 
modulates  the  flux  of  this  energy  to  the  deeper  water  masses.  Conversely, 
flux  of  heat  back  to  the  atmospheric  boundary  layer  has  an  important 
influence  upon  the  climate  and  its  fluctuations.  _ 

Garwood  (1976)  and  Zilitinkevich  et  al.  (1979)  provide  background  summaries  of  historical 

contributions  to  the  field  of  mixed-layer  dynamics.  Early  modeling  frameworks  involved 

only  the  mean  momentum  budget,  with  no  inclusion  of  mechanical  energy  -  either  buoyancy 

or  turbulent  kinetic  energy.  Kraus  and  Turner  (1967)  first  included  mechanical  stirring  in 

their  model,  parameterizing  shear  production  of  turbulent  kinetic  energy  (TKE)  on  the  wind 

stress.  Niiler  (1974)  emphasized  the  importance  of  using  the  turbulent  kinetic  energy 

equation  in  modeling  the  mixed  layer,  and  included  entrainment  shear  production  and 

damping  in  the  TKE  budget.  Garwood  (1977)  refined  the  mechanical  energy  budgets  used 

so  that  entrainment  depended  on  the  vertical  energy  component  rather  than  bulk  total 

turbulent  energy.  In  a  more  detailed  three-component  TKE  equation  model,  mixed-layer 

deepening  was  shown  to  depend  on  wind  direction,  Garwood  et  al.  (1985). 

Price  et  al.  (1986)  used  a  simplified  one-dimensional  bulk  gradient  Richardson 

number-based  model  to  study  the  separate  and  combined  effects  of  diurnal  heating  and 

cooling,  wind  stress,  and  entrainment.  The  model  gave  realistic  profiles  of  temperature  and 
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velocity,  and  predicted  a  trapping  depth  for  the  upper  ocean’s  response  to  the  diurnal  forcing 
cycle. 

McDougall  (1987)  coined  the  term  “thermobaricity”  in  his  study  of  the  effects  of 
nonlinearities  in  the  seawater  equation  of  state  for  isopycnal  mixing.  Gill  (1973)  first 
included  thermobaricity  in  a  static  model  for  shelf  convection.  The  dynamic  mixed-layer 
importance  of  thermobaricity,  which  is  the  pressure  dependence  of  the  thermal  expansion 
coefficient,  which  he  termed  thermobaricity,  was  explained  by  Garwood  (1991),  using  a 
model  that  individually  calculated  the  three  components  of  turbulent  kinetic  energy. 

Large  et  al.  (1994)  have  studied  the  use  of  a  scheme  called  K  Profile  Parameterization 
(KPP)  which  uses  a  nonlocal  transport  term  K  for  the  boundary  layer,  assigns  an  eddy 
viscosity  K  for  the  ocean  interior,  and  matches  the  K  profile  at  the  boundary  layer  interface. 
This  approach  does  allow  a  representation  of  entrainment,  but  the  nonlocal  transport  term  is 
determined  only  by  surface  fluxes  and  does  not  parameterize  the  effect  of  entrainment  fluxes. 

2.  Deep  Convection 

Early  observations  in  the  Mediterranean  Sea  show  that  regions  of  deep  water 
formation  are  selected  by  the  large-scale  stratification,  to  first  order  (Gascard,  1973).  Later 
Killworth  ( 1 979, 1 983)  noted  that  chimneys,  or  regions  of  open-ocean  deep  water  formation, 
are  narrower  than  can  be  explained  by  large-scale  gyres  and  cooling  regimes;  some  other 
preconditioning  mechanism  must  determine  the  specific  area  for  overturning.  Cyclonic 
eddies  produced  by  baroclinic  instability  of  the  mean  flow  tend  to  reduce  stratification  by 
lifting  and  cooling  a  localized  area  of  approximately  chimney  scale  (~20km);  such  eddies 
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lifting  and  cooling  a  localized  area  of  approximately  chimney  scale  (~20km);  such  eddies 
were  Killworth’s  likely  suspects  for  the  preconditioning  mechanism. 

From  observations  in  the  Labrador  Sea,  Clarke  and  Gascard  (1983)  found  open-ocean 
deep  water  formation  occurring  under  approximately  the  same  conditions  specified  by 
Killworth.  Further  observations  made  in  winter  1994-1995,  reported  in  Lilly  et  al.  (1998), 
showed  intermediate-depth  winter  convection,  with  200-1000  m  wide  plumes  reaching 
1750m  in  depth.  The  depth  varied  greatly  with  location,  and  the  water-mass  properties  were 
not  horizontally  homogenized  to  that  depth. 

Greenland  Sea  observations  in  the  winter  of  1988-1989  highlight  the  variable  nature 
of  convection  in  the  Greenland  Sea.  There  is  seasonal  variability,  as  the  area  becomes  ice- 
covered  and  then  ice-free  again  when  the  Nordbukta  opens,  and  interannual  variability,  as 
winter  forcing  produces  mixing  to  great  depth  in  some  years,  and  only  intermediate  depths 
in  others  (Schott  et  al.,  1993).  In  1989,  Morawitz  et  al.  (1996)  measured  a  chimney  of  scale 
40  km  or  less,  which  penetrated  to  below  1000m,  and  speculated  that  its  location  was 
preconditioned  by  a  hole  in  the  underlying  Arctic  Intermediate  Water,  resulting  from 
unsteadiness  of  inflow  in  the  summer  of  1 988.  However,  a  net  salinity  gain  for  the  column 
was  observed  but  unexplained.  Observations  in  the  winter  of  1993-1994  also  showed  a 
feature  that  could  be  classified  as  a  chimney,  an  anticyclonic  eddy  with  horizontal  scale  ~50 
km.  Although  convection  was  limited  to  intermediate  depth  throughout  the  basin,  this 
feature  was  the  site  for  the  deepest  mixing  observed  (800m).  However,  the  winter  of  these 
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observations  was  an  atypical  one,  since  no  Is  Odden  formed  and  the  surface  buoyancy  loss 
was  less  than  the  decadal  mean  (Lherminier,  1998). 

Marshall  and  Schott  (1999)  provide  a  recent  summary  of  observations,  theory  and 
models  of  open-ocean  deep  convection,  including  discussions  of  the  importance  of  rotation, 
thermobaricity,  and  horizontal  inhomogeneity.  They  note  that  the  interplay  of  plume  scales, 
mesoscale  eddy  scales,  and  large-scale  preconditioning  is  important  in  all  known  sites  of 
open-ocean  deep  convection. 

3.  Various  Modeling  Approaches 

The  importance  of  deep  convection  to  the  climate  has  motivated  many  modeling 
approaches.  Ultimately,  climate  models  must  be  able  to  realistically  convey  heat  and  tracers 
from  the  atmosphere  to  the  intermediate  and  bottom  waters  of  the  ocean.  An  important  step 
in  this  complex  problem  is  the  development  of  a  computationally  simple  one-dimensional 
model  that  can  be  applied  to  each  vertical  column  of  a  GCM  grid,  and  result  in  realistic 
vertical  profiles  and  water  mass  production. 

The  simplest  representation  of  the  polar  atmosphere’s  effect  on  the  heat  distribution 
of  the  ocean  is  that  buoyancy  loss  at  the  surface  is  compensated  by  overturning  to  remove 
hydrostatic  instability  through  convective  adjustment.  Sander  et  al.  (1995)  compared  the 
results  of  a  nonhydrostatic  model  to  those  for  pure  convective  adjustment.  He  concluded  that 
convective  adjustment  may  be  adequate  when  small-scale  dynamics  are  relatively 
unimportant,  but  that  the  vertical  advective  transport  of  tracers  is  inhibited.  In  such  schemes. 
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buoyancy  loss  is  conceived  as  a  mechanism  for  overturning  by  convective  adjustment  only, 
and  is  not  recognized  as  an  important  source  of  mechanical  energy  for  entrainment. 

Visbeck  et  al.  (1995)  used  an  ice-coupled  model  to  study  the  factors  necessary  for 
deep  water  formation  in  the  Greenland  Sea.  They  found  that  a  strictly  one-dimensional 
approach  could  not  account  for  the  differences  in  deepening  between  the  water  columns  of 
die  Nordbukta  and  those  of  the  Is  Odden.  However,  a  quasi  two-dimensional  approach  that 
included  an  average  ice  export  did  accurately  reproduce  observed  ocean  response. 

Garwood  et  al.  (1994)  hypothesized  two  kinds  of  thermobaric  conditional 
instabilities:  a  thermobaric  layer  instability,  in  which  a  layer  may  become  statically  unstable 
if  advected  downward;  and  a  thermobaric  parcel  instability,  in  which  parcels  overshoot  the 
mixed-layer  interface,  become  thermobarically  unstable,  are  detrained  from  the  mixed  layer, 
and  accelerated  downward.  Detraining  plumes  were  simulated  in  an  LES  model,  and 
proposed  scalings  were  given  for  critical  depth  and  velocity.  Garwood  and  Isakari  (1993) 
developed  a  scheme  for  use  in  ocean  general  circulation  models,  in  place  of  the  convective 
adjustment  routine,  that  includes  the  effects  of  parcel  instabilities,  with  possible  detrainment 
of  fluid  from  the  mixed  layer.  Paluskiewicz  and  Romea  (1997)  report  on  the  development 
of  a  similar  scheme. 

Hakkinen  (1992)  used  a  large-scale  three-dimensional  coupled  ice-ocean  model  with 
2.5  order  turbulence  closure  (Mellor  and  Yamada,  1974),  without  a  convective  adjustment 
scheme,  to  model  deep  convection  in  the  Greenland  Sea.  Hakkinen  imposed  a  barotropic 
mean  flow  over  shelf  break  topography  containing  either  a  canyon  or  a  spur.  She  found  that 
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in  the  presence  of  an  ice  edge,  which  induces  upwelling,  either  the  canyon  or  the  spur  can 
provide  preconditioning  and  determine  a  deep  convection  site. 

Several  studies  focus  not  on  the  initiation  and  location  of  deep  convection,  but  on  its 
later  development  and  effects  once  a  given  region  has  been  subject  to  deep  mixing  (Jones  and 
Marshall,  1993;  Legg  and  Marshall,  1993;  Send  and  Marshall,  1995;  Jones  and  Marshall, 
1997).  These  studies  hypothesize  that  the  sinking,  spreading,  and  restratification  phases  of 
deep  convection  will  be  modeled  effectively  by  largely  neglecting  plume-scale  processes, 
representing  them  by  a  vertical  mixing  time  scale,  and  focusing  on  the  breakup  of  the  mixing 
site  by  baroclinic  instability. 

Stone  (1997)  used  a  one-dimensional,  three-component  model  to  study  mixed-layer 
deepening  in  the  Greenland  Sea  as  observed  during  February  and  March  1 994.  Comparison 
of  energy  contributions  from  buoyancy  flux  versus  wind  forcing  showed  a  surprising 
dominance  of  wind  energy,  and  underscored  the  need  for  further  study  of  the  relative 
importance  of  wind  and  cooling  in  models  of  deep  mixed  regimes. 

Akitomo  (1999a,  1999b)  investigated  the  effects  of  thermobaricity  on  open-ocean 
deep  convection.  High-resolution  numerical  experiments  with  constant  eddy  viscosity  and 
eddy  diffusivity  were  initialized  with  realistic  ocean  conditions;  rapid  overturning  occured 
in  the  Weddell  Sea  due  to  thermobaric  instability,  and  more  gradual  deepening  occured  in 
the  Greenland  Sea  due  to  greater  stratification. 

Yoshikawa  et  al.  (1999)  used  a  large-scale  nonhydrostatic  model  with  constant  eddy 
viscosity  to  study  convective  processes  in  the  presence  of  a  geostrophic  velocity  shear.  They 
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found  that  the  convective  plumes  accelerate  the  growth  of  baroclinc  instability  and  enhance 
vertical  motion. 

Straneo  and  Kawase  (1999)  compared  two  models  of  localized  convection,  one  due 
to  localized  forcing,  and  the  other  due  to  uniform  forcing  but  localized  preconditioning. 
They  found  that  the  localized  forcing  tends  to  overemphasize  the  importance  of  baroclinic 
instability  during  active  convection,  and  increasing  the  horizontal  gradients  during  the  course 
of  convection,  whereas  convection  in  the  preconditioned  scenario  decreased  the  horizontal 
gradients  initially  present.  Studies  using  a  two-dimensional  model  (Straneo  et  al.,  1999) 
show  that  the  interaction  between  horizontal  gradients  and  active  convection  results  in 
slantwise  convection  and  can  alter  mixed-layer  deepening  from  that  predicted  by  one¬ 
dimensional  considerations.  Marshall  and  Schott  (1999)  provide  a  discussion  of  slantwise 
convection,  and  discuss  the  transition  from  convection  to  baroclinic  instability  for  cases  with 
localized  forcing. 

The  focus  of  this  dissertation  is  a  better  understanding  of  the  physics  involved  in  deep 
mixing;  there  are  still  gaps  in  the  detailed  knowledge  of  ocean  turbulence  that  may  affect 
large-scale  parameterizations.  An  important  limitation  to  many  models  being  used  to  study 
deep  convection  is  the  use  of  hydrostatic  physics  and/or  the  use  of  unrealistic  eddy 
viscosities  and  eddy  diffusivities  even  when  nonhydrostatic  terms  are  resolved.  Detailed 
data  sets  of  geophysical  turbulence  are  needed  to  test  how  well  viscous  models  apply  to  the 
ocean.  Because  of  the  difficulty  of  obtaining  data  sets  of  geophysical  turbulence,  large  eddy 
simulation  is  an  important  tool  for  improving  our  understanding  of  deep  turbulent  mixing 


11 


in  the  ocean.  We  turn  to  large  eddy  simulation  because  of  its  ability  to  represent  turbulence 
to  as  fine  a  scale  as  desired;  the  grid  spacing  is  chosen  small  enough  to  resolve  energy  down 
to  scales  including  a  part  of  the  -5/3  log-slope  inertial  energy  cascade.  Processes  on  scales 
smaller  than  this  are  parameterized  (see  Chapter  II,  Section  A). 

4.  Large  Eddy  Simulation 

Large  eddy  simulation  (LES)  is  a  modeling  technique  that  resolves  the  energy- 
containing  and  flux-supporting  eddies  at  and  below  the  integral  scales  of  motion.  Unlike 
direct  numerical  simulation  (DNS),  LES  does  not  resolve  the  Kolmogoroff  microscale,  and 
eddies  smaller  than  the  LES  resolution  are  parameterized.  Unlike  Reynolds  average 
numerical  simulation  (RANS),  the  parameterization  uses  an  eddy  viscosity  that  is  time  and 
space  dependent  (Orszag,  1993).  Smagorinsky  (1963)  and  Lilly  (1967)  pioneered  the  use 
of  such  nonlinear  eddy  viscosities  to  parameterize  the  effect  of  unresolved  turbulence. 
Deardorff  (1973)  added  a  dynamic  equation  for  the  subgrid  energy  to  determine  the  nonlinear 
eddy  viscosity.  In  an  atmospheric  boundary  layer  model,  he  was  able  to  achieve  realistic 
temperature  variances,  turbulent  intensities,  and  to  preserve  the  sharpness  of  the  inversion 
base  during  entrainment.  Moeng  (1984)  created  a  new  LES  model  for  study  of  the 
atmospheric  planetary  boundary  layer,  following  Deardorff.  She  used  spectral  methods  in 
the  horizontal  to  take  advantage  of  FFT  algorithms  and  improved  closure  assumptions  in  the 
subgrid  parameterization  by  using  a  spectral  cut-off  filter  to.separate  large-eddy  and  subgrid 
domains.  Further  work  by  Moeng  and  Wyngaard  (1988)  refined  the  subgrid  dissipation  and 
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eddy  viscosity  parameterizations  to  improve  the  inertial  range  spectra.  They  also 
demonstrated  the  relative  insensitivity  of  the  model  to  the  subgrid  parameterization. 

The  atmospheric  boundary  layer  LES  model  of  Moeng  (1984)  was  modified  to 
represent  the  oceanic  boundary  layer,  and  was  used  to  simulate  escaping  thermobaric  plumes 
by  Garwood  et  al.  (1994).  The  work  included  free  convection  in  a  neutral  layer,  and  free 
convection  in  a  neutral  layer  with  salinity-stratified  water  below.  Some  plumes  that  reached 
the  interface  passed  through  and  were  termed  thermobaric  plumes  because  the  increase  in  the 
thermal  expansion  coefficient  caused  them  to  change  from  being  positively  buoyant  at  the 
interface  to  neutrally  and  then  negatively  buoyant  as  the  parcels’  momentum  carried  them 
below  the  interface. 

The  LES  model  has  been  further  modified  to  improve  the  resolved  spectra  by  using 
horizontally  isotropic  filtering  methods  (Harcourt  and  Garwood,  1994),  and  to  control  errors 
by  using  an  upwind  numerical  scheme  for  scalar  advection  near  sharp  boundaries  (Harcourt, 
1999). 
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II.  THEORETICAL  OVERVIEW 


The  ocean  mixed  layer,  characterized  both  by  turbulence  and  by  temperatures  and 
salinities  relatively  homogenized  with  depth,  resides  above  the  nonturbulent,  stably  stratified 
interior  of  the  ocean.  A  sharp  density  gradient  is  usually  at  the  lower  boundary  of  the  mixed 
layer,  while  its  upper  boundary  is  the  wavy  air-sea  interface.  For  the  highly  variable  fluxes 
from  the  atmosphere  to  have  any  effect  upon  the  underlying  water  masses,  they  must  be 
communicated  through  this  turbulent  layer. 

Although  temperature  and  salinity  are  characteristically  unstratified  within  the  mixed 
layer,  the  amount  of  inhomogeneity  varies  greatly.  Variances  within  the  mixed  layer  depend 
on  surface  fluxes  of  temperature  and  salinity,  via  cooling,  rain,  and  entrainment  of 
contrasting  fluid  from  below.  Lateral  mixing  in  a  horizontally  inhomogeneous  regime  also 
produces  variance. 

The  physical  properties  of  seawater  itself  add  complexity  to  mixed-layer  dynamics. 
Near  freezing  temperature,  the  thermal  expansion  coefficient  is  especially  sensitive  to 
pressure;  this  effect,  termed  thermobaricity,  can  produce  instabilities  when  a  cold  fresh 
mixed  layer  overlies  warmer,  more  saline  water.  This  type  of  instability  may  result  in 
detrainment  of  mixed-layer  water  into  the  stratified  water  below.  Another  interesting 
consequence  of  thermobaricity  is  the  possibility  of  fluids  of  equal  density  but  contrasting 
temperature  and  salinity  laterally  interleaving,  with  the  warm,  salty  intrusions  being 
preferentially  preserved. 
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A.  FUNDAMENTAL  EQUATIONS 


1.  General  Equations 

The  Navier-Stokes  equations  of  motion  are  the  starting  point  for  calculating  mixed- 
layer  turbulence: 
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The  indices  i  =  1,2,3  correspond  to  the  Cartesian  coordinate  directions  of  east,  north,  and 
up;  Q  is  earth’s  rotation:  [Q.]  =.  Q[  0  cos<(>  sin<j)  ];  molecular  viscosity  is  v; 
gravitational  acceleration  is  g ;  and  <}>  is  latitude. 

The  turbulent  kinetic  energy  (TKE)  equation  is  formed  from  the  equations  of  motion 

by  decomposing  the  variables  into  their  mean  and  turbulent  parts,  ut  =  u.  +  «/ ,  where  the 


turbulent  part  is  not  assumed  to  be  small  compared  to  the  mean.  Using  the  Boussinesq 
approximation  and  incompressibility,  the  total  TKE  equation  for  the  system  is 
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The  equation  for  mean  flow  is  formed  from  equation  (2.A.1)  by  decomposing, 
averaging,  making  the  Boussinesq  approximation,  and  neglecting  viscous  effects  on  the  mean 


flow: 
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The  equations  for  scalars  potential  temperature  and  salinity  are  formed  from 


conservation  of  energy  and  mass: 
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where  kt  and  ks  are  molecular  diffusivities. 

2.  Turbulent  Energy  Cascade 

It  is  important  in  a  turbulence-resolving  model  to  be  sure  that  the  smallest  scales 
resolved  are  part  of  the  inertial  subrange.  In  an  energy  spectrum,  if  the  energy  per  unit 
wavenumber  E  (m3  s'2)  is  a  function  of  only  the  wave-number  k  (m'1)  and  viscous  dissipation 

e  (m2  s*3),  then  the  energy  in  the  larger  scales  cascades  down  to  smaller  scales  following 
E  «  c2/3  k  5/3  (2.A.6) 

The  portion  of  the  spectrum  that  follows  this  -5/3  law  is  referred  to  as  the  turbulent  inertial 
subrange,  and  connects  the  scales  at  which  the  system  is  being  forced  to  the  scales  at  which 
dissipation  is  occurring  (Kolmogoroff,  1941). 

3.  One-Dimensional  Equations 

A  set  of  one-dimensional  equations  derived  from  the  more  general  equations  (2.A.1)- 
(2.A.5)  is  presented  briefly  here;  full  details  of  its  derivation  can  be  found  in  Stone  (1997). 
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The  general  equations  have  been  decomposed,  averaged,  and  vertically  integrated  over  a 
depth  h.  The  hydrostatic  and  Boussinesq  approximations  have  been  made;  horizontal 
homogeneity  has  been  invoked;  and  viscous  effects  on  the  mean  flow  have  been  neglected. 


The  equations  for  the  ageostrophic  mean  flow  are 


=  2Q3h<V> 
dt  3 


(2.A.7) 


a(/?<K>)  =  -  2 Q,h<U>  +  ^  (2.A.8) 

dt  p0 

The  geostrophic  balance  has  been  subtracted  from  these  horizontal  momentum  equations. 
Overbars  indicate  averaging  in  time,  the  angle  brackets  indicate  the  vertical  average,  and  the 
wind  stresses  are  t  and  t  . 

x  y 

The  temperature  and  salinity  equations  are 

4^  =  i(  A.  -  w  A0  )  (2.A.9) 

dt  h  P0cp 

T-  =  T<  *Fsr  ~  we )  (2-A-10) 

at  h 

The  temperature  flux  at  the  surface  has  been  written  as 

0V|  0  =  —  (2. A.  11) 

PocP 

with  Q0  the  rate  of  cooling  of  the  surface,  and  cp  the  specific  heat  of  seawater  at  constant 
pressure.  The  surface  salinity  flux  is 

(2.  a.12) 

with  Fsw  being  the  net  rate  of  equivalent  water  fluxed  upward  at  the  surface  due  to 
precipitation,  evaporation,  melting,  and  freezing,  and  having  units  meters/second.  Fluxes  of 
conserved  quantities  at  the  bottom  of  the  mixed  layer  are  calculated,  assuming 
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approximately-linear  flux  profiles  within  the  layer,  as 


C'w' |  =  -w  AC 


(2.A.13) 


where  C  is  the  conserved  quantity  and  we  is  the  entrainment  velocity.  Also,  —  =  we 
in  the  absence  of  upwelling  or  downwelling. 


The  turbulent  kinetic  energy  equations  are 
dh<Un>  _  +  +  jQiji 


dt 


u*  P0 


_____  I 

+  2  mJ  <E>  -  3<«  /2>  )<E>2  -  —m<E> 

3 


(2.  A.  14) 


dh<v  > 

It 


2  m,  T  „  _ 

- (— )  +  w  e(Av)2 


W 


Po 


+  2m,(  <E>  -  3 <v /2>  )<£>2  -  —m.<E> 

3 


(2.A.15) 


dh<w/2>  ,  r  fio  ,  ai^  .  R  —  T,  ,  , 

- - -  =  gh[  - (  a0+  -  )  +  p  S  Vol  ] 

dt  pc  3 


-  ghwe[  A0(  a0+  -2—  )  -  P  AS  ]  -  2 Q2h  — 

1-5  P0 

_  j_  2 

+  2m J  <E>  -  3<w/2>  )<E>2  -  - m.<E> 2  (2.A.16) 

3 

The  thermal  expansion  coefficient  is  treated  not  as  constant,  but  is  allowed  to  vary  linearly 
with  depth: 

a(z)  =  a0  -  QjZ  (2A.17) 

The  constants  mvm2,m3  are  the  dissipation  constant,  the  pressure  redistribution  constant, 
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and  the  shear  production  constant.  The  friction  velocity  is  given  by  u  *2  =  —  J  t2  +  t2  , 

_  _  _ Po V 

the  total  turbulent  kinetic  energy  is  <E>  =  (  <u/2>  +  <v/2>  +  <w 72 >  ),  and  the 

3_ 

dissipation  is  D  =  2mi<E>2  . 

4.  Scalings 

Some  useful  scalings  are  presented  here  for  use  in  later  chapters.  First,  a  measure  of 
buoyant  production  of  TKE  due  to  surface  fluxes  is  the  free-convection  velocity  scale  w0* 
calculated  from  surface  fluxes  in  the  TKE  equation 

(2.A.18) 

Jo  h  pc  3 

r  p 

This  is  the  rate  at  which  TKE  is  produced  by  hydrostatic  adjustment  of  buoyancy  loss.  If  ctj 
is  zero,  then 

v,;3  =  -B0h  (2.A.19) 

When  a,  is  nonzero,  the  velocity  scale  includes  the  additional  energy  produced  by  downward 
acceleration  of  parcels  due  to  thermobaricity.  The  total  buoyant  production  is  given  by  a 
velocity  scale  that  includes  the  effect  of  buoyant  damping 

w*3  =  -2fh^Y'dz  (2.A.20) 

Jo 

Since  w0*3  is  known,  and  w*3  can  be  calculated,  the  buoyant  damping  of  TKE  due  to 
entrainment  can  be  deduced 


*3  *3 

WE  =  w 


(2.A.21) 


A  natural  Rossby  number  for  a  layer  of  depth  h  in  free  convection,  found  by  scaling  of  the 


momentum  equation,  is 

R  =  —  (2.A.22) 

°  fh 

The  Rossby  internal  radius  of  deformation,  the  distance  a  parcel  must  travel  before  the  effect 
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of  rotation  becomes  significant,  is 

_  {7jL 


Rd  =  (2.A.23) 

and  can  be  applied  to  a  layer  homogeneous  in  potential  temperature  but  with  depth- 
dependent  a(z)  as 

(2.A.24) 


*<  = 


f 


B.  EQUATIONS  FOR  SCALAR  VARIANCE 


Beginning  with  the  temperature  equation  (2.A.4),  and  the  equation  for  mean 
temperature  formed  by  decomposition  and  averaging  equation  (2.A.4), 

(2.B.1) 


39  =  -,,'221  ♦  vV>5 

dt  '  dxt  T 

equation  (2.B.1)  is  subtracted  from  equation  (2.A.4)  to  give  an  equation  for  the  temperature 


■30' 


perturbation. 


/30' 


30/  30  /u\J  Til  a! 

-  =  -u. -  +  It, -  +  K-V2©' 

dt  'a*.  dx.  T 


Multiplying  equation  (2.B.2)  by  0'  and  averaging  gives 

=  -  ^ 
dt  2  '  dx. 


u.  0'—  -  m/0'—  +  Kr  0'V2©' 


dx. 


dx. 


(2.B.2) 


(2.B.3) 


If  we  assume  that  0/2  is  horizontally  homogeneous,  w  =  0  ,  and  neglect  curvature 
effects  on  diffusion,  equation  (2.B.3)  becomes 


1/2 


d,Q'\  ~^de  td  / 0/2  3  ,0'2  3  ,0\ 

—  ( - )  =  -  u  0 -  -  w  0 -  -  ( - U  -  +  - V  -  +  — w  - ) 

dt  2  dx  dz  dx  2  dy  2  dz  2 


K  30^  30^  +  30^30^  +  30'  30' 

r  dx  dx  dy  dy  dz  dz 


(2.B.4) 
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The  first  two  terms  on  the  right-hand-side  are  horizontal  and  vertical  gradient  production,  the 
first  term  in  parentheses  is  turbulent  transport,  and  the  last  is  dissipation. 

For  application  to  the  one-dimensional  model  of  this  paper,  we  require  horizontal 
homogeneity  of  0  and  integrate  Equation  2.B.4  over  a  mixed  layer  of  depth  h  .  Net 
turbulent  transport  vanishes  since  there  is  no  flux  of  variance  through  the  surface  or  the 
interface.  Since  the  vertical  derivative  of  0  is  non-zero  only  near  the  surface  and  near  the 
interface,  the  integral  of  the  gradient  production  term  can  be  written  as  two  separate  terms. 
The  contribution  from  the  interface  is 


We(A0)2 


(2.B.5) 


where  w  is  the  entrainment  velocity  and  A0  is  the  potential  temperature  jump  across 


the  interface.  The  contribution  from  the  surface  is 

-  60  =  (—)(——)  =  — (— f 


(2.B.6) 


P CP  u*  pcp 


w  Pc„ 


Q 

where  — —  is  the  surface  heat  flux,  m3  is  a  constant  of  proportionality,  and  the  near- 

PS 

surface  gradient  60  depends  on  both  the  surface  heat  flux  and  a  velocity  scale  w*  in  forced 
convection,  or  w*  in  free  convection.  The  dissipation  rate  of  temperature  variance  must 
depend  on  the  amount  of  variance  present  and  on  some  dissipation  time  scale  ,  such  that 


X  H  Kr 


a07  607 

dxj  dxj 


m. 


<0/2> 

T 

x 


(2.B.7) 


where  m,  is  a  constant  of  proportionality,  and  angle  brackets  mean  the  vertical  average  has 
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been  taken.  Variance  must  dissipate  on  the  same  time  scale  on  which  energy  dissipates,  that 
h 

is,  x  =  —  ,  so 


f°  X  dz  =  mj<0/2>  <E> 

J  -h 


1/2 


The  one-dimensional  temperature  variance  budget  then  has  the  form 

d_ 
dt 


ik^j  -  -  .,<¥>  <*>« 


u*  Pc 


(2.B.8) 


(2.B.9) 


In  steady  state,  solving  for  temperature  variance  gives  the  result 


w3(  Qo  )2  +  W,(A0)2 


_  u  mc.  2 


<Q'Z>  = 


Pc„ 


mi<E> 


1/2 


(2.B.10) 


The  first  term  in  the  numerator  is  the  surface  production  of  temperature  variance,  and  the 
second  is  entrainment  production. 

Similarly,  if  the  salinity  flux  at  the  surface  is  negligible, 


<S,2>  =  e 


w  (AS)2 


1/2 


2m1<E> 


C.  THERMOBARIC  INSTABILITIES 


(2.B.11) 


The  thermal  expansion  coefficient  of  seawater, 


a 


=  _  1  dP 

Po  00 


varies  with  both  temperature  and  pressure. 


(2.C.1) 


The  variation  with  pressure  can  give  rise  to 


instabilities  that  are  termed  thermobaric  instabilities.  The  vertical  dependence  of  a  can  be 
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expressed  as 

da 

a(z)  =  a0-  ?-z 
dz 


(2.C.2) 


The  thermobaric  depth  scale  Ha  indicates  the  depth  at  which  the  thermal  expansion 


coefficient  is  twice  its  surface  value 
-cco 

H  =  — — 
a  8a 

dz 


dec 

Taking  —  as  a  constant,  a(z )  can  be  written 
dz 

a(z)  =  a0(l  - 

Hcc 

1.  Thermobaric  Plumes 


(2.C.3) 


(2.C.4) 


A  parcel  that  is  displaced  downward,  such  as  in  Figure  (2),  will  experience  a  change 
in  buoyancy 

A  b  =  -ga1d1Az  (2.C.5) 

doL 

where  z.  =  —  .  If  this  buoyancy  difference  is  greater  than  the  environment’s  buoyancy 
dz 

difference  for  the  two  levels,  the  parcel  will  continue  its  vertical  motion  away  from  its 
original  position.  Mixed  layer  parcels  near  the  interface  experience  displacements  downward 
during  active  convection,  and  if  the  displacement  az  is  greater  than  a  critical  depth 


h  =  (  ft---  -1  )Ha  (2.C.6) 

cr  aoA0 

thermobaric  plumes  may  escape  from  the  mixed  layer  (Garwood  et  al.,  1994).  The  critical 


depth  would  be  reached  by  a  parcel  with  initial  downward  speed 

a0£A6  1/7 

Wcr=(hcr~  > 


(2.C.7) 
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2.  Thermobaric  Interleaving 

It  is  possible  for  a  parcel  to  be  laterally  displaced  in  such  a  way  that  it  forms  a 
temperature  and  salinity  contrast  with  the  parcels  immediately  above  and  below  it,  but  leaves 
the  profile  of  potential  density  unperturbed,  such  as  in  Figure  (3).  The  density-compensated 
intrusion  has  temperature  and  salinity  characteristics  such  that 

<V0&6  =  P  &S  (2.C.8) 


and  the  background  buoyancy  gradient  is 

db  ^0  n  dS  ,To  .  «  .  /A  y~. 

—  =  ga—  ~  g P—  =  N2  (z  0  )  (2.C.9) 

dz  dz  dz 

We  define  the  vertical  coordinate  r|  =  0  at  the  level  of  the  intrusion,  and 

« =  v„o  -■?->  C-C.io) 

If  the  intrusion  is  displaced  vertically,  it  experiences  a  vertical  acceleration  due  to  buoyant 


forces 

=  b'  =  agQ'  - 
dt2 

where 

07  =  60  -  ii  — 
dz 

and 


(2.C.11) 


(2.C.12) 


S'  =  5S  -  r\—  (2.C.13) 

dz 

Substituting  equations  (2.C.8),  (2.C.9),  (2.C.10),  (2.C.12),  and  (2.C.13)  into  (2.C.1 1),  we  get 


the  wave  equation 


dt 2 


g  a6/ 

(- — +  iV2)q  =  0 
H „ 


(2.C.14) 


with  stable,  oscillating  solutions  for  warm  intrusions  (0'  >  0).  For  cold  intrusions 


(0X  <  0)  there  is  a  critical  background  stratification 
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below  which  the  intrusion  will  be  unstable  to  vertical  perturbations.  Thus  we  expect  that  in 
regions  where  thermobaricity  is  strong,  profiles  may  reveal  some  preference  for  warm  salty 
intrusions,  and  that  cold  fresh  intrusions  will  be  less  commonly  observed. 


minutes. 


III.  MODELS  FOR  LARGE  EDDY  AND  ONE-DIMENSIONAL  SIMULATION 


A.  ONE-DIMENSIONAL  FORMULATION 

The  one-dimensional  equations  of  Chapter  II  are  closed  by  scaling  the  TKE  budget 

at  the  bottom  of  the  mixed  layer.  Shear  production  and  dissipation  are  small  compared  to 

the  transport  of  turbulence  from  above  and  the  buoyant  damping  due  to  entrainment,  so  that 

the  balance  is  mainly  among  storage,  turbulent  transport,  and  buoyant  damping.  See  Stone 

(1997)  for  details.  The  entrainment  velocity  we  is 
m .  <E>  V  <w,2> 

",  =  -r-i - - - —  (3.A.1) 

<E>  +  g/i(aA0  -  PAS') 

The  dimensionless  constant  m4  is  introduced  to  compensate  for  the  fact  that  one- 
dimensional  equations  cannot  account  for  the  vertical  distribution  of  energy  within  the  layer, 
and  turbulence  is  typically  more  concentrated  in  the  upper  mixed  layer  than  near  the 
entrainment  zone. 

Table  1  summarizes  the  one-dimensional  model  equations  used  in  this  study.  The 
system  of  equations  was  cast  as  a  MATLAB  function,  and  the  function  was  then  called  using 
the  MATLAB  utility  ODE45,  a  Runge-Kutta  ordinary  differential  equation  solver  that  has 
automatic  step-size  control. 


27 


ae  =  j_(  Qq 

St  h  p  ecf 


-  w#A0  ) 


—  =  i(  SFsw  -  W'AS  ) 
dt  h  sw  ’ 

-2Q3h<v>  +  ^ 
St  3  P„ 


*<*<F>>  =  -  20,  h<v>  +  A 

St  ^  p0 

m4  <£>  \/  <w/2> 

W  -  — = - — - — 

<E>  +  gh( aA0  -  PAS) 


+  w,(A«)2  +  20^/1  — 
dt  «,  p0  *  P0 


2m2(  <£>  -  3<u /2>  )<E>2  -  jm1<£>2 


=  — (^-)2  +  *,(Av)2 
5/-  W.  p0 


7  2 


2w2(  <£>  -  3<v  >  )<E> 2  -  “Wj<£> 


Bh<w,2> 


ar 


gA[  -£L(  ao+  fill  )  +  p  5  Vol  ]  +  2m2(  <£>  -  3<w,2>  )<£>2  -  -ra,<£>2 

3  3 

-  ghwt[  A0(  a0+  — -j-  )  -  (3  AS  ]  -  2 CJjA  — 


Table  1.  One-Dimensional  Model  Equations. 
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B.  LARGE  EDDY  FORMULATION 


The  LES  model  calculates  three-dimensional,  nonhydrostatic,  geophysical  turbulence 
using  the  vorticity  forms  of  equation  (2.A.1)  with  the  Boussinesq  approximation,  the 
continuity  equation,  plus  the  heat  and  salinity  budgets  expanded  from  equation  (2.A.4): 

Bu  Bp  *  dx  3t  3t  (3.B.T) 

dt  dx  y  dx  dy  dz 


Bv  Bp  *  St  dx  dx 

—  =wC  -«C  -2Q  u— SL—2L—Z. 

dt  dy  z  dx  dy  dz 


(3.B.2) 


Bw  Bp  *  St  dx  dx 


(3.B.3) 


du  +  dv  +  dw  _ 
dx  dy  dz 


(3.B.4) 


S6  30  30  30  3  30,  A  3  30.  3  30, 

dy  dz  dx  dx  dy  6  dy  dz  0  dz 


dt 

dS 


dx 


dx 
dS . 


dS  dS  dS  d  35,  3  35,  3  35, 

dy  dz  dx  s  dx  dy  s  dy  dz  s  dz 


(3.B.5) 


(3.B.6) 


dt  dx  dy  dz  dx  *  dx' 

Here  u,  v,  and  w  are  the  easterly,  northerly,  and  vertical  velocity  components,  C  is  vorticity,  2  Q 


is  the  vertical  Coriolis  parameter,  2  Qy  is  the  horizontal  Coriolis  parameter,  and  a  is  the 

thermal  expansion  coefficient.  Shear  stresses  are  expressed  as  xy  such  that 
du.  du. 


(3.B.8) 


The  dynamic  pressure  P*  includes  terms  for  resolved  and  unresolved  energy  e 


Po  3 

In  the  model  code,  a  is  calculated  using  a0  and  a,  as  in  equation  (2. A.  1 7).  Specified  values 
for  wind  and  surface  buoyancy  flux  enter  the  equations  through  the  surface  boundary 


00 

conditions  by  specifying  values  for  xi3  in  the  momentum  equations  and  values  for  K%  —  and 

dz 

pj  rr 

Ks —  in  the  scalar  equations.  The  bottom  boundary  condition  is  slip  with  respect  to  the 
dz 

mean  flow  and  no-slip  with  respect  to  perturbations  from  the  mean.  Boundary  conditions  are 
doubly  periodic  in  x  and  y.  The  system  of  equations  is  solved  using  second  order  centered 
finite  differencing  in  the  vertical,  the  spectral  method  of  Fox  and  Orszag  (1973)  in  the 
horizontal,  and  time  advancement  with  an  Adams-Bashforth  scheme. 

Turbulence  is  assumed  to  be  isotropic  at  the  grid  scale  and  smaller.  Subgrid  scale 
fluxes  of  momentum,  salinity,  and  temperature  are  calculated  using  second  order  turbulence 
closure,  giving  eddy  mixing  coefficients  KM ,  Ks ,  and  KQ ,  respectively,  that  are  time-  and 
space-dependent  (Smagorinsky,  1963).  The  coefficients  depend  on  the  length  scale  X  of  the 
unresolved  turbulence,  the  grid  scale  6  and  the  velocity  scale  given  by  the  unresolved 
turbulent  kinetic  energy  e: 

km  =  Ke  =  Ks  =  ck  (3-B-9> 

The  subgrid  turbulence  length  scale  in  neutral  and  unstable  conditions  is  the  resolution  scale 

6,  related  to  the  grid  spacing  Ajc,  Ay,  and  A z  by: 

X  =  8  =  (c6  Ax  Ay  A z)m  (3.B.10) 
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The  equation  for  unresolved  TKE  includes  terms  for  advection,  shear  production,  buoyancy 


flux,  turbulent  transport,  and  dissipation  (Deardorff,  1973): 

£  -  -  *  tK4 3  -  p  *)  ♦  -  « 

‘dx.  "  dx.  5  6  dz  v  dz  dx.  Mdx  ’ 

where 

du  du 

x.:  =  -  i  +  --*■) 

"  dx/ 

Subgrid  dissipation  £  is  modeled  as: 

•  e3/2 

€  =  C, - 


(3.B.11) 


(3.B.12) 


(3.B.13) 


When  stratification  is  stable,  the  length  scale  of  the  unresolved  turbulence  may  be 
smaller  than  the  resolution  scale  6,  so 

X  =  min(6,Xstab)  (3.B.14) 

where 

fe 


^ stab  ^ stab 


N 


(3.B.15) 


and  N  is  the  buoyancy  frequency.  The  dissipation  constant  is  also  reduced  in  stable 
conditions,  using: 


c€  =  +  ce2g 


(3.B.16) 

Choice  of  model  constants  cK  =  0. 128,  c6  =  1 . 14,  cstab  =  0.76,  ce  =  2.29,  c6,  =  0. 1 9,  and 
=2.1  follows  Harcourt  (1999).  Additional  details  concerning  the  basic  numerical  method 
are  provided  by  Moeng  (1984),  Garwood  et  al.  (1994),  and  Harcourt  (1999). 
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C.  IMPOSITION  OF  HORIZONTAL  GRADIENTS  IN  THE  LES 


Since  it  was  desirable  for  the  LES  model  to  retain  doubly  periodic  boundary 
conditions,  zonal  gradients  were  imposed  in  the  LES  domain  not  by  changing  temperature 
or  salinity  across  the  box  explicitly,  but  by  adding  a  gradient  term  to  the  x  derivative  after  it 
is  calculated  in  the  spectral  domain  and  passed  back  to  the  spatial  domain.  The  zonal 
advection  terms  of  equations  (3.B.5)  and  (3.B.6)  are  calculated 

—  =  ...(os  previously  calculated). ..+  tgrad  (3.C.1) 

dx 

=  ...(as  previously  calculated).. .  +  sgrad  (3.C.2) 

dx 

The  pressure  gradient  term  for  u  momentum,  equation  (3.B.1),  has  an  additional  term  that 


is  depth  dependent 

=  ...(as  previously  calculated)...*  pgrad(iz)  (3.C.3) 

dx 

The  depth  dependent  pressure  gradient  is  calculated  starting  at  the  bottom  level  in  the  LES 


domain  and  integrating  upward  to  produce  the  correct  thermal  wind  shear 

pgrad(nnz)  =  -CL  tgrad  —  +  P  sgrad  — 

2  2 

pgrad(iz)  =  pgrad(iz  +  1)  -  a  tgrad  dz  +  (5  sgrad  dz 


(3.C.4) 

(3.C.5) 


where  a  and  P  are  computed  level  by  level. 
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IV.  GREENLAND  SEA  ENTRAINING  MIXED  LAYER  CASE  STUDY 


Previous  work  recommended  the  comparison  of  the  one-dimensional  model  and  LES 
for  Greenland  Sea  entrainment  (Stone,  1997).  Balance  among  energy  components  and  the 
dominance  of  wind  or  cooling  are  particularly  at  issue,  and  the  one-dimensional  model 
constants  have  been  adjusted  to  better  emulate  LES-predicted  balances. 

A.  DESCRIPTION  OF  PARAMETERS 

Both  the  LES  model  and  the  one-dimensional  model  were  initialized  with 
temperature  and  salinity  data  from  a  sounding  taken  by  the  R/V  Valdivia  near  74.5°N,  2.5°W 
on  February  16th,  1994.  There  was  a  200-m  deep  mixed  layer  of  cold  fresh  water  over 
relatively  warm  and  saline  water.  Another  sounding  was  taken  on  March  19th,  31  days 
later;  the  mixed  layer  had  deepened  to  600  m.  The  profiles  used  as  the  initial  conditions 
were  subjected  to  a  convective  adjustment  scheme  to  avoid  initial  instabilities  below  the  200- 
m  layer.  Comparison  of  heat  and  salinity  content  of  the  water  column  between  the  two 
soundings  showed  an  average  loss  of  130  W/m2  and  no  net  salinity  flux.  Wind  and  heat  flux 
data  from  the  European  Center  for  Medium  Range  Weather  Forecasts  (ECMWF)  were  used 
to  force  both  models,  reduced  by  a  factor  of  0.56  so  that  the  average  heat  loss  over  the  period 
matched  that  observed  in  the  water  column.  Applicability  of  this  one-dimensional 
assumption  was  discussed  in  Stone  (1997).  Figure  4  shows  the  wind  speed  and  net  heat 
fluxes  for  the  period  simulated,  as  well  as  the  predicted  mixed-layer  depths.  Values  for  % 
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and  a,  were  3.6e-5  and  2.8e-8,  respectively. 

Values  for  the  dissipation,  pressure  redistribution,  shear  production,  and  entrainment 
constants  used  in  the  one-dimensional  model  are  shown  in  Table  2. 


m. 

m2 

m* 

m4 

1.0 

3.0 

1.9 

0.3 

Table  2.  Values  of  One-Dimensional  Model  Constants 

The  LES  simulation  was  done  on  a  grid  having  50  vertical  levels,  96x96  horizontal 
grid  points,  20m  horizontal  grid  spacing  and  1 9m  vertical  grid  spacing.  The  grid  spacing  for 
horizontal  and  vertical  grids  differs  so  that  the  resulting  resolution  is  isotropic.  The  strong 
wind  forcing  led  to  high  surface  velocities  and  required  a  short  time  step  of  30  seconds. 

B.  COMPARISON 


The  comparison  of  results  from  one-dimensional  and  LES  models  focuses  on  the 
relative  balance  of  energy  components  and  the  relative  importance  of  TKE  budget  terms. 
Figure  6  shows  TKE  components  from  each  model.  The  total  TKE  is  similar,  and  although 
both  models  show  that  TKE  is  anisotropic,  the  energy  is  distributed  differently  among 
components  for  the  two  models.  The  one-dimensional  model  has  most  energy  in  the  v 
component  during  the  peaks  because  most  of  the  strong  wind  forcing  was  from  southerly 
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winds.  The  LES  energy  components  tend  to  distribute  wind  shear  production  more  evenly 
between  u  and  v  because  it  resolves  shear  and  turning  between  each  vertical  grid  level. 


Figures  7  and  8  show  u  and  v  TKE  budget  terms;  there  are  significant  differences  in 
the  sources  of  turbulent  energy.  The  one-dimensional  model  equations  include  the  three 
TKE  component  equations,  so  the  terms  of  the  TKE  budget  plotted  in  Figures  7  and  8  are 
explicit  within  the  one-dimensional  model.  The  LES  TKE  production  terms  plotted  are 


du 


dv 


calculated  from  the  mean  covariances  u  w  and  v  w  and  the  mean  shears  —  and  —  for 


dz  dz 

shear  production.  In  general,  the  one-dimensional  model  is  sensitive  to  wind  direction 
through  the  shear  production  term  while  the  LES  allows  shear  to  spiral  as  it  is  transmitted 
downward,  so  that  shear  production  from  v  wind  is  put  into  both  u  and  v  turbulent  energy. 
The  pressure  redistribution  term  of  the  one-dimensional  model  is  the  term  that  moves  energy 
produced  by  v  wind  shear  into  u  and  w  components,  but  the  resulting  balance  is  not  the  same 
as  that  in  the  LES. 

Figure  9  allows  comparison  of  buoyant  production  and  buoyant  damping.  Buoyant 
production  due  to  surface  forcing  for  both  models  is  computed  as  w0*2  using  Equation 
2.A.19.  Buoyant  damping  for  the  one-dimensional  model  is  calculated  using  the  term  in 
Equation  2.A.16  that  contains  we.  To  calculate  LES  total  buoyancy  flux  w*3,  either  the  mean 
covariances  w  'Q1  and  w 's or  the  net  mean  changes  in  0  and  S  profiles  over  many  timesteps 
can  be  used.  In  this  case  the  change  in  mean  profile  method  was  chosen.  Buoyant  damping 
for  the  LES  is  then  computed  as  total  buoyancy  flux  minus  buoyant  production  due  to  surface 
flux,  as  shown  in  Equation  2. A.2 1 . 
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The  ratio  of  entrainment  buoyancy  flux  to  surface  buoyancy  flux  used  by  Stull  (1988) 


■Ar 


-b  'w  'I 


-h 


b  'w  '\ 


is  applied  here  to  the  ocean  with  thermobaricity  in  a  more  general  form, 

*3 


/  _ 


(4.B.1) 


(4.B.2) 


Note  that  A  '  R  =  AR  only  for  an  idealized  bulk  mixed  layer  in  the  absence  of  thermobaricity. 
The  ratios  are  averaged  over  the  entire  period  and  over  the  two  periods  with  the  strongest 
cooling  and  are  shown  in  Table  3.  These  values  are  large  compared  to  standard  values  for 
free  convection  (Stull,  1988)  due  both  to  the  strong  wind  stress  and  to  thermobaricity 
(Harcourt,  1999). 


Entire  Period 

(days  47.7-83.8) 

Days  53.0-57.5 

Days  74.0-78.0 

1-D 

.42 

.35 

.44 

LES 

.36 

.40 

.31 

Table  3.  Comparison  of  ratios  of  entrainment  damping  to  buoyant  production. 


The  one-dimensional  model  plot  indicates  that  during  the  strongest  forcing  event,  occurring 
on  about  Julian  day  76,  vertical  TKE  is  getting  about  as  much  energy  from  wind  shear  via 
pressure  redistribution  as  from  cooling. 

Figure  10  shows  the  total  TKE  budget  terms  for  each  model.  Shear  production  is  not 
as  dominant  in  the  LES  as  in  the  one-dimensional  model,  but  is  still  2-3  times  as  large  as 
buoyant  production.  These  balances  are  among  bulk  quantities,  however;  shear  production 
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in  the  one-dimensional  model  is  immediately  available  for  entrainment,  while  in  the  LES 
much  of  the  energy  produced  by  shear  is  produced  in  the  upper  levels  and  dissipated  locally. 
Figure  11  is  similar  to  Figure  10  except  the  buoyant  production  and  damping  terms  are 
summed  to  show  the  net  effect. 
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V.  SIMULATIONS  OF  THERMOBARIC  PARCEL  INSTABILITIES 


The  thermobaric  plumes  described  in  section  2.C.1  were  simulated  using  the  LES 
model  for  mixed-layer  and  interior  stratification  conditions  designed  to  be  similar  to 
wintertime  Greenland  Sea  stratification  conditions.  Parcels  first  detrained  from  the  mixed 
layer  in  agreement  with  the  predictions  for  critical  depth  hcr  and  critical  velocity  wcr 
hypothesized  by  Garwood  (1994)  in  Equations  2.C.6  and  2.C.7.  The  escaping  plumes  exhibit 
a  characteristic  mushroom-cap  shape,  and  affect  water  properties  at  depth. 

A.  DESCRIPTION  OF  PARAMETERS 


The  simulations  of  parcel  instabilites  began  with  a  180m-deep  mixed  layer  over  a 
layer  of  neutral  or  stably  stratified  fluid.  Uniform  cooling  of 200  W/m2  and  5  m/s  northward 
wind  was  applied  at  the  surface.  At  the  onset  of  cooling,  the  sharp  interface  was  perturbed 
with  inertial  waves  excited  by  convective  plumes,  and  entrainment  began.  All  the 
simulations  used  a  grid  spacing  of  10m,  128x128  horizontal  grid  points,  and  80  vertical 
levels,  except  the  case  with  the  strongest  stratification,  which  had  the  lowest  30  levels 
removed  to  save  computing  time. 

Two  groups  of  experiments  were  conducted;  in  the  first  group,  rotation  and 
thermobaricity  were  varied.  Table  4  shows  experiment  designations  and  parameters.  The 
case  titled  Rotation  and  Full  Thermobaricity  included  Qy  and  Qz  and  calculated  depth- 
dependent  a  throughout.  The  case  titled  No  Rotation  had  Qy  =  Qz  =  0.  The  case  titled 
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Thermobaricity  Off  Below  the  Mixed  Layer  included  rotation  terms,  but  calculated  a  as 
depth-dependent  down  to  220m,  then  kept  the  220m  value  for  the  fluid  below.  In  this  way, 
turbulence  within  the  mixed  layer  retained  the  effect  of  thermobaricity,  but  parcels  at  the 
interface  were  not  accelerated  downward  if  they  reached  depths  below  220m.  Initial  mixed- 
layer  temperature  and  salinity  properties  were  -1 . 13  °C  and  34.834  psu;  lower  layer  properties 
were  -0.43  °C  and  34.900  psu.  Values  for  a0  and  ax  were  6.5e-5  and  2.75e-8,  respectively, 
giving  an  initial  buoyancy  drop  at  the  base  of  the  mixed  layer  of  3.4142e-6  m/s2. 


Rotation  and  Full 
Thermobaricity 

No  Rotation 

Thermobaricity  Off 
Below  the  Mixed  Layer 

gsea3 

gsea4 

gsea5 

Table  4.  First  group  of  thermobaric  plume  simulations,  all  with  unstratified  lower  layer. 


Unstratified 

N2  =  0 

Weak 

Stratification 

N2  =  8.5e-8 

Moderate 

Stratification 

N2  =  1.7e-7 

Strong 

Stratification 

N2  =  3.9e-7 

gsea3.2 

gstrat2 

gstrat4 

gstrat 

Table  5.  Second  group  of  thermobaric  plume  simulations,  with  varying  stratification  below 
mixed  layer. 


In  the  second  group  of  experiments,  stratification  below  the  mixed  layer  was  varied. 
Table  5  shows  experiment  designations  and  parameters.  In  each  case,  the  mixed-layer 
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temperature,  salinity,  and  the  buoyancy  jump  at  the  base  of  the  mixed  layer  were  intialized 
identically.  Figure  12  shows  the  initial  temperature  and  salinity  profiles;  Figure  13  shows 
the  stratification  in  the  various  experiments.  Figure  14  shows  the  relative  buoyancy  a 
mixed-layer  parcel  would  have  if  plunged  to  any  depth  below  the  initial  interface.  The  zero 
intercept  of  each  relative  buoyancy  profile  is  the  depth  below  which  a  mixed-layer  parcel 
would  become  hydrostatically  unstable.  Values  for  a0  and  a,  were  3.69e-5  and  2.65e-8, 
respectively,  and  the  initial  buoyancy  jump  varied  only  slightly  from  that  in  the  first  group 
of  experiments,  at  3.4289e-6  m/s2.  Representative  vertical  velocity  spectra  are  shown  in 
Figures  15  and  16. 

B.  PLUME  FORMATION 

Detrainment  of  parcels  occurred  in  all  the  unstratified  cases  and  in  weak  and 
moderate  stratification.  Detection  of  escaping  plumes  was  accomplished  using  graphical 
animation.  Figure  17  shows  a  typical  field  of  branching  and  mushroom-cap-shaped 
detraining  plumes. 

Figure  18,  top  frame,  contains  a  plot  of  mixed-layer  depth  for  the  unstratified 
experiments.  Mixed-layer  depth  was  determined  by  the  maximum  gradient  of  the 
horizontally-averaged  temperature,  spline  interpolated  between  grid  levels.  The  critical 
depth  hcr  of  Equation  2.C.6  is  also  plotted.  Since  the  interface  is  vigorously  deformed  and 
oscillating,  the  bottom  of  the  mixed-layer  should  be  considered  the  median  position  of  a 
sharp  but  moving  interface,  rather  than  as  the  location  of  strongest  gradient  within  a  diffuse 
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region.  The  mixed  layer  is  deepening  and  becoming  more  saline  through  entrainment  and 
is  cooling  due  to  the  surface  heat  loss.  The  simulation  without  rotation  entrains  slightly 
faster,  without  the  stabilizing  effect  of  rotation.  Critical  depth  is  reduced  because  the  ratio 
is  decreasing.  The  buoyancy  jump  across  the  interface,  which  was  initially  salinity- 

a0A6 

dominated,  is  becoming  less  so  as  it  weakens.  After  detrainment  begins,  the  mixed  layer 
shallows  as  a  result  of  the  mass  lost  from  the  mixed  layer  by  detrainment.  The  standard  case 
was  run  to  catastrophic  detrainment,  that  is,  until  escaping  fluid  volume  had  completely 
depleted  the  upper  layer. 

The  vertical  velocities  plotted  in  Figure  18  are  the  rms  values  for  the  bulk 
horizontally-averaged  w  variance.  The  vertical  velocity  variance  is  in  steady  state,  but  the 
critical  velocity  wcr  of  Equation  2.C.7  is  decreasing.  Plume  formation  occurs  for  both  of  the 
thermobaric  cases  when  wcr  drops  to  within  about  twice  the  rms  fluctuation  of  w. 

The  bottom  frame  in  Figure  18  shows  the  buoyant  production  averaged  over  the 
domain  below  the  mixed  layer,  which  is  normally  zero  without  plume  detrainment.  Buoyant 
production  is  calculated  from  the  covariance  of  vertical  velocity  perturbations  with 
temperature  and  salinity  perturbations,  plus  subgrid  contributions.  The  experiment  with 
rotation  off  has  plumes  escaping  earliest,  since  vertical  motions  are  not  constrained.  Of  the 
two  predictors  for  the  onset  of  escaping  plumes,  hcr  and  wcr,  the  critical  depth  exhibits  more 
clearly  the  difference  in  timing  or  the  start  of  detrainment  between  the  experiments  with  and 
without  rotation;  each  begins  buoyant  production  below  when  the  mean  position  of  the 
mixed-layer  interface  is  about  50m  above  the  critical  depth. 
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Figure  19  shows  4-hour-averaged  profiles  of  buoyant  production  before  and  during 
detainment,  except  for  the  case  with  thermobaricity  off,  which  had  no  detainment,  but  is 
plotted  for  comparison.  Buoyant  production  within  the  mixed  layer,  which  normally  has  a 
nearly-linear  profile,  is  affected  by  the  buoyant  flux  profile  below. 

Figure  20,  top  frame,  shows  a  plot  of  mixed-layer  depth  and  for  the  set  of  stratified 
experiments.  For  these  experiments,  the  mixed  layer  also  deepens  and  becomes  more  saline 
through  entrainment  and  cools  due  to  surface  cooling.  Critical  depth  shallows  as  the 
buoyancy  jump  across  the  interface  weakens.  The  moderate  stratification  case  was  continued 
until  catastrophic  detainment  occured. 

The  center  frame  in  Figure  20  depicts  the  rms  values  for  the  bulk  horizontally- 
averaged  w  variance  and  the  critical  velocities.  The  vertical  velocity  variance  is  near  steady 
state,  but  wcr  is  again  decreasing  as  h,  hcr,  and  A0  all  change.  Plume  detainment  does  occur 
when  wcr  is  less  than  twice  the  rms  fluctuation  of  w,  but  buoyant  production  below  the  mixed 
layer,  shown  in  the  bottom  frame  of  Figure  20,  is  not  as  clear  an  indicator  for  escaping 
plumes.  This  is  because  there  is  strong  buoyant  damping  occurring  as  the  plumes  entrain 
water  from  their  new  surrounding  stratified  environment. 

Figure  21  shows  profiles  of  buoyant  production  before  and  during  detainment. 
Buoyant  production  within  the  mixed  layer  is  relatively  unaffected  by  the  process  below. 
Comparing  Figures  19  and  21,  buoyant  production  again  is  seen  to  be  a  less  clear  indicator 
of  detainment  in  the  stratified  cases  than  it  is  in  neutral  stratification.  Figure  22  plots  the 
skewness  of  vertical  velocity 


The  magnitude  of  this  parameter  is  small  for  wave  motion  (  SKw « 1  ).  Increases  in 
magnitude  of  SKW  indicate  that  large  deviations  from  the  mean  are  concentrated  on  one  side 
of  the  mean;  the  distribution  is  skewed.  In  the  LES-produced  data,  skewness  at  a  given  level 
changes  suddenly  with  the  passage  of  plumes  through  that  level.  A  depth  of  285m  was 
chosen  for  evaluation  because  it  was  not  affected  by  the  mixed-layer  interface,  and  because 
there  was  no  significant  delay  between  plume  detrainment  and  the  arrival  of  plumes  at  that 
depth.  A  skewness  of  -0.7  was  chosen  somewhat  arbitrarily  for  use  in  the  following  section 
as  a  consistent  indicator  for  the  onset  of  detrainment. 


C.  PLUME  ENTRAINMENT  AND  EFFECT  ON  INTERIOR  FLUID 


The  effect  of  detraining  plumes  on  the  interior  of  the  ocean  can  be  seen  in  Figure  23 
for  the  unstratified  set  of  experiments.  With  thermobaricity  included  at  all  depths,  there  is 


more  effect  on  temperature  in  the  lower  levels  than  the  mid-levels.  For  the  experiments 
where  stratification  is  varied,  the  effect  is  not  illustrated  well  by  the  color  plots  of 
temperature  (Figure  24).  Profiles  of  temperature  before  and  during  detrainment  are  similarly 
difficult  to  compare  (Figure  25).  Since  there  is  no  salinity  flux,  the  fraction  of  mixed-layer 
water  present  at  any  level  can  be  calculated  using  the  initial  salinity  at  that  level,  the  current 


salinity  at  that  level,  and  the  current  mixed-layer  salinity: 

S(z )  ,  -  S(z). , 

A  (z)  ~  V  Jcurrent  K  'initial  (5  C  1) 

Figure  26  plots  the  fraction  of  mixed-layer  water  present  at  each  level,  10  and  20  hours  after 


detrainment  begins.  This  quantity  better  distinguishes  the  differences  among  the  scenarios. 
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In  the  environments  with  nonzero  stratification,  the  effect  of  plume  entrainment  did  not 
produce  a  region  at  depth  that  contained  a  greater  fraction  of  mixed-layer  water  than  the 
fraction  present  at  mid-depths  below  the  mixed  layer. 

The  parameterization  of  Paluskiewicz  and  Romea  (1997)  assumes  a  top-hat 
distribution  for  a  plume’s  temperature  and  salinity  across  a  horizontal  slice,  and  determines 
the  plume  entrainment  rate  as  proportional  to  a  calculated  vertical  velocity  profile,  which 
typically  has  a  maximum  below  the  mixed  layer.  While  no  direct  comparison  of 
detrainment’s  effect  on  interior  water  properties  as  calculated  by  the  parameterization  and 
by  LES  is  made  here,  a  check  against  LES  data  for  these  two  assumptions  made  by  the 
scheme  can  be  done.  Figures  27-30  show  temperature  and  salinity  cross-sections  for 
detraining  plumes.  Temperature  and  salinity  are  relatively  constant  across  the  plumes  in  all 
simulations,  with  gradients  concentrated  at  a  sharp  edge,  but  more  so  for  the  unstratified 
cases.  Figures  31  and  32  shows  vertical  cross-sections  of  w  that  capture  several  detraining 
plumes.  Plume  velocities  can  have  local  maxima  below  the  mixed  layer  in  both  the  stratified 
and  unstratified  cases. 
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VL  SIMULATIONS  IN  THE  PRESENCE  OF  A  LARGE-SCALE  GRADIENT 


Active  convection  and  postconvective  spin-down  were  modeled  with  horizontal 
gradients  imposed  in  the  east-west  direction  using  the  scheme  of  Section  3.C.  The  horizontal 
gradients  result  in  mean  north-south  currents  generated  by  thermal  wind  shear.  Mean 
profiles  of  u  show  a  slow  overturning  trend.  Scalar  variances  of  temperature  and  salinity  and 
TKE  increase  as  a  result  of  the  large-scale  gradient.  Furthermore,  horizontal  fluxes  and 
stability  within  the  mixed  layer  are  affected  by  the  interaction  of  overturning  and  turbulence. 


A.  DESCRIPTION  OF  PARAMETERS 


The  numerical  experiments  in  this  chapter  used  a  domain  2.56km  x  2.56km  x  400m 
deep,  128x128  horizontal  grid  points,  and  20  vertical  levels,  with  the  exception  of  the  case 
with  no  surface-forced  convection.  This  case  was  1 .28km  x  1.28km  x  400m  deep  and  had 
41  vertical  levels.  Aside  from  the  mean  horizontal  gradients,  the  LES  domain  was  initially 
horizontally  and  vertically  homogeneous,  so  that  when  convection  began,  the  entire  domain 
became  fully  turbulent,  with  no  entraining  layer.  Uniform  cooling  of  either  200  or  400  W/m2 
was  applied,  with  negligible  surface  wind  stress.  Table  6  shows  the  parameters  that  were 
varied  in  the  different  experiments.  The  latitude  of  75°N  determined  Qz ,  but  Qy  was  held  to 
a  zero  value.  The  suffix  “tbo”  on  an  experiment  designation  indicates  that  the  effect  of 
thermobaricity  was  turned  off  by  setting  a}  to  zero. 
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Numerical  Experiments  -  LES  and  Horizontal  Stratification 

Surface  Forcing  Regime  || 

Horizontal 

Gradient 

Regime 

Quiescence 

Active 

Convection 

Qo=200  W/m2 

Active 

Convection 

Q0=400  W/m2 

Post 

Convection 

Qo=0 

No  Gradient 

gr5200 

gr5,  gr5tbo 

gr5co 

Weak  Gradient 

gradbcomp 

gr9 

gr4,  gr4tbo 

gb7,  gb7tbo 

Strong  Gradient 

gr6 

gr6co 

Table  6.  Names  and  titles  of  horizontal  gradient  experiments. 


The  two  nonzero  horizontal  gradient  regimes  are  described  by  the  parameters  in  Table 
7,  and  Figure  31  plots  the  buoyancy  gradients  and  resultant  geostrophic  velocities.  In  both 
regimes,  temperature  and  salinity  decrease  eastward,  and  their  gradients  are  partially 


6 

B 

a® 

a, 

tgrad 

sgrad 

max 

bgrad 

Rd 

°C 

■■ 

°C'1 

°C-'  m'1 

'Cm'1 

psu  m'1 

s’2 

m 

Weak 

Gradient 

-1.5 

34.85 

3.15e-5 

2.92e-8 

-6.25e-6 

-6.05e-7 

2.7e-9 

1835 

Strong 

Gradient 

-1.0 

34.85 

3.84e-5 

2.86e-8 

-1.33e-4 

-6.70e-6 

-1.2e-8 

1573 

Table  7.  Values  used  in  the  horizontal  gradient  experiments. 
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compensating  in  buoyancy.  In  the  weak  gradient  experiment,  salinity  (sgrad)  dominates  the 
buoyancy  gradient  ( bgrad ),  and  a  northward  geostrophic  flow  results.  The  buoyancy  gradient 
changes  with  depth  due  to  the  pressure  dependence  of  the  thermal  expansion  coefficient,  in 
conjunction  with  the  temperature  gradient  ( tgrad ). 

bgrad  =  ~g(aQ  +  a  xz)  tgrad  +  g  ft  sgrad  (6.A.1) 

The  small  buoyancy  gradient,  with  temperature  and  salinity  gradients  almost  completely 
compensating,  was  originally  chosen  in  hopes  of  generating  the  intrusions  discussed  in 
Section  2.C.2.  No  cases,  either  with  or  without  cooling,  exhibited  interleaving. 


*d  = 


The  vertical  density  gradient  is  used  to  compute  the  Rossby  radius  of  deformation 
=  1 -  (6.A.2) 


/ 


The  sizes  of  the  deformation  radii  are  close  to  the  box  dimension  of  2.56km;  features  the 
scale  of  the  Rossby  radius  cannot  be  well  resolved  in  these  numerical  experiments. 

In  the  strong  gradient  experiment,  properties  and  gradients  were  chosen  to  be  similar 
to  a  rim  current  region  at  the  edge  of  the  Greenland  Sea  gyre  near  75°N  and  6°E  in  1989, 
reported  by  Budeus  et  al.  (1993),  with  the  thermal  expansion  coefficient  chosen  to 
correspond  to  a  temperature  of  1.0  °C.  Temperature  and  salinity  decrease  eastward,  but  the 
effect  of  temperature  dominates  the  buoyancy  gradient.  This  buoyancy  gradient  is  an  order 
of  magnitude  larger  than  in  the  weak  gradient  experiment.  Thermal  wind  shear  produces 
surface  currents  on  the  order  of  1-2  cm/s  for  both  cases,  but  in  opposite  directions.  Steady 
state  was  reached  in  the  cases  with  no  gradient  and  weak  gradient,  but  a  growing  horizontal 
box  mode  was  present  in  the  strong  gradient  case;  spectra  are  shown  in  Figures  32  and  33. 
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B.  SIMULATIONS  OF  ACTIVE  CONVECTION 


A  particularly  striking  feature  of  the  simulations  with  large-scale  gradients  is  die 
appearance  of  features  such  as  those  shown  in  the  salinity  field  of  Figures  34.  The  horizontal 
variability  is  organized  at  scales  larger  than  the  scale  of  convective  plumes,  shown  in  Figure 
35.  Temperature  variability,  shown  in  Figure  36,  is  similarly  organized  into  larger  scales, 
but  the  appearance  of  convergence  zones  of  sinking  surface-cooled  water  partially  masks  the 
organization  into  larger  scales  in  the  near-surface  temperature  field. 

1.  Bulk  Turbulent  Kinetic  Energy  and  Variance 

Figure  37  shows  the  bulk  mean  TKE  for  experiments  with  varying  gradients.  The 
total  amount  of  TKE  increases  with  the  strength  of  the  horizontal  gradient.  Most  of  the 
increase  is  in  the  horizontal  TKE  components;  Figure  38  shows  the  ratio  of  vertical  to  total 
TKE;  the  ratio  decreases  with  the  strength  of  the  horizontal  gradient.  Figure  39  shows  bulk 
temperature  and  salinity  variances.  For  the  different  experiments,  the  amount  of  variance 
increases  with  the  strength  of  the  horizontal  gradient.  The  increase  in  TKE  increases  subgrid 
diffusivity,  but  the  continual  advection  of  more  disparate  parcels  into  the  domain  results  in 
a  net  increase  in  the  variance.  Experiments  with  weak  and  no  large-scale  gradients  showed 
little  change  in  time,  but  in  the  strong  gradient  experiment,  the  variance  increase  was  evident 
even  before  the  TKE  showed  significant  increase.  The  strong  gradient  experiment  did  not 
achieve  steady  state  even  after  a  period  of  more  than  16  days.  The  peaks  in  TKE  and 
variance  were  associated  with  the  growth  and  destruction  of  a  feature  at  the  scale  of  the  LES 
domain,  shown  in  Figures  40-42.  The  nature  of  this  feature  was  quite  complicated:  clear 
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signatures  were  found  in  the  temperature,  salinity,  buoyancy,  and  pressure  fields,  but  not  in 
the  vertical  velocity  field;  and  the  horizontal  velocity  field  shows  detailed  structure. 

Table  8  lists  some  common  quantities  used  in  scaling  turbulence  and  compares  the 
results  in  these  experiments  with  the  scalings  suggested  by  Harcourt  (1999).  These  scalings 

Vii\ 

suggest  that - =  0.375  i?02495  and - s  0.558  R0Z9S  .  Compared  with  these 

w  *2  w  *2 

scalings,  the  values  found  here  are  consistently  lower  by  about  ten  percent,  showing 
disagreement  in  the  scaling  coefficient.  Although  w  *  increases  slightly  with  larger  horizontal 
gradients,  the  scaling  relation  between  w /2  and  w*2  appears  to  be  less  sensitive  to  the 
horizontal  gradients.  The  scaling  relations  were  determined  by  Harcourt  (1999)  using 
simulations  without  the  thermobaric  effect;  the  results  in  Table  8  show  slightly  better 
agreement  in  the  cases  with  thermobaricity  off. 


Natural  Rn 

w'2/TKE1/2 

IffTHSl 

B9 

gr5 

2.34 

.44 

.42 

.29 

.31 

.39 

.44 

gr5tbo 

2.34 

.44 

.42 

.28 

.31 

.40 

.44 

gr5200 

1.86 

.35 

.42 

.27 

.29 

.36 

.41 

gr4 

2.34 

.44 

.40 

.29 

.31 

.39 

.44 

gr4tbo 

2.38 

.44 

.40 

.29 

.31 

.41 

.44 

gr9 

1.86 

.35 

.40 

.27 

.29 

.37 

.41 

mm 

2.48 

.46 

.33 

.29 

.31 

.41 

.45 

Table  8.  Turbulent  scaling  quantities  for  the  experiments  with  horizontal  gradients. 
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2. 


Currents 


Figure  43  shows  the  time-mean  profiles  of  the  horizontal  velocity  components.  The 
north-south  velocity  is  similar  to  the  geostrophic  balance  required  for  the  buoyancy  gradient 
that  was  prescribed.  Velocities  vanish  at  the  bottom  of  the  domain  because  the  frame  of 
reference  for  the  LES  moves  with  the  mean  flow  at  the  bottom  grid  level.  The  mean  east- 
west  current  profile  for  the  weak  gradient  cases  shows  westward  flow  near  the  surface,  which 
would  bring  in  fresher  and  cooler  water  above,  and  eastward  flow  near  the  bottom,  which 
would  advect  in  warmer  and  more  saline  water  below.  Since  the  horizontal  buoyancy 
gradient  is  salinity-dominated,  the  surface  inflow  is  relatively  buoyant,  and  the  lower-level 
inflow  is  more  dense.  This  tendency  for  more  buoyant  water  to  be  advected  in  near  the 
surface  and  less  buoyant  water  to  be  advected  in  near  the  bottom  is  a  feature  of  overturning. 
Overturning  currents  move  the  system  toward  a  more  stable  stratification  state.  In  the  strong 
gradient  case,  the  mean  east- west  current  profile  shows  eastward  flow  in  most  of  the  layer, 
which  would  advect  in  warmer  and  more  saline  water  above,  and  westward  flow  near  the 
bottom,  which  would  advect  in  cooler  and  fresher  water  below.  The  buoyancy  gradient  is 
temperature-dominated,  so  again  the  system  is  overturning.  Comparison  of  the  overturning 
currents  for  the  cases  with  weak  gradient  and  varying  amounts  of  cooling  shows  that 
convection  increases  the  speed  with  which  the  system  advects  fluid  in  the  effort  to  reach  a 
state  with  lower  potential  energy.  In  a  rotating  reference  frame  such  as  this,  mean  u  currents 
affect  mean  v,  so  resolved  values  for  v  are  different  from  the  geostrophic  values.  Figure  44 
shows  schematically  this  effect  at  the  first  grid  level. 
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Figure  45  shows  time  plots  of  the  mean  velocity  at  the  first  grid  level.  Pure  inertial 

oscillation  for  the  latitude  75°N  has  a  period  of  12.5  hours,  and  the  case  with  no  convection 

exhibits  this  period  of  oscillation.  The  cases  with  cooling  show  that  convection  does  have 

an  effect  on  the  horizontally-averaged  inertial  currents.  The  period  is  modulated:  the  weaker 

cooling  case  shows  oscillations  with  a  period  of  1 1  hours;  the  strong  cooling  cases  show 

oscillations  with  periods  of  10-11  hours.  Also,  the  coherence  of  the  oscillation  in  the 

horizontal  mean  flow  at  the  surface  varies  in  time,  first  being  disrupted,  but  then 

reestablishing  a  coherent  oscillation  as  ageostrophic  motions  continue  to  perturb  the  steady- 

state  balance.  These  effects  were  not  investigated  in  detail,  but  the  time  scales  of  the 

turbulent  eddies  were  close  to  the  inertial  period,  so  that  interaction  between  convective  and 

H 

inertial  motions  seems  likely.  Convective  time  scales  -  were  calculated  as  1 1 .6  hours 

Wm,s 

for  the  weak  gradient  -  200  W/m2  case,  9.0  hours  for  the  weak  gradient  -  400  W/m2  case,  and 
8.3  hours  for  the  strong  gradient  -  400  W/m2  case. 

3.  Horizontal  Fluxes 

Mean  currents  in  a  horizontally  nonhomogeneous  regime  advect  salinity  and 
temperature.  Because  mean  currents  are  really  the  result  of  time-averaging  a  turbulent  flow, 
the  advection  by  the  mean  flow  can  be  resolved  as  turbulent  fluxes  in  the  LES.  Figure  46 
shows  the  mean  profiles  of  horizontal  fluxes  for  the  weak  and  no  gradient  cases. 
Considering  fluxes  in  the  direction  of  the  geostrophic  mean  flow  first,  v  'Q'  and  v's' ,  there 
is  net  movement  of  cold  temperatures  northward  in  the  weak  gradient  cases,  as  the  mean 
northward  current  carries  away  newly  cooled  water.  Net  movement  of  more  saline  water 
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northward  is  intensified  near  the  bottom  because  the  overturning  continually  moves  more 
saline  water  into  the  lower  region. 

There  are  also  turbulent  fluxes  caused  by  the  overturning.  Salinity  flux  u  's 1  is  bow¬ 
shaped  as  westward  current  near  the  surface  brings  fresher  water  in  with  u  /0/  >0  ,  and 
eastward  current  near  the  bottom  brings  in  more  saline  water,  again  with  u  /0/  >0  . 
However,  u ' 0 '  is  mostly  negative  for  the  case  with  weak  gradients  and  400  W/m2  cooling. 
This  means  there  is  net  westward  transport  of  warm  temperatures  despite  the  overturning 
current’s  tendency  to  move  cooler  water  west  in  the  upper  part  of  the  column.  In  the 
experiment  with  weak  gradients  and  200  W/m2  cooling,  the  turbulent  flux  u  fr7  is  close  to 
zero  over  most  of  the  column,  where  again  a  positive  value  would  be  expected.  The  reasons 
for  this  unexpected  result  are  unclear,  but  could  be  due  in  part  to  the  fact  that  the  turbulent 
part  of  velocity,  u ' ,  is  computed  by  subtracting  the  mean  for  each  level,  not  the  mean  over 
the  entire  domain,  or  due  in  part  to  the  effective  shifting  of  the  domain  as  the  mean  flow  at 
the  lowest  level  is  subtracted.  In  any  case,  the  effect  of  this  flux  in  the  variance  budget  of 
Section  2.B  turns  out  to  be  relatively  small.  Using  Equations  2.B.4  and  2.B.10,  the  steady- 
state  one-dimensional  bulk  variance  budget  predicts  that  the  temperature  variance  is 

-2  fh  ^Q'^-dz  -  2  f  * 

Jo  dz  Jo  dx 

<6>  =  - 

mx<E>m 

=  2 _ ^ _ A1  -  2  f  Vo'dz  (6.B.1) 

w*mx<E>V2  Pcp  m1<E>m  j° 

The  first  term  on  the  right-hand-side  is  vertical  gradient  production,  and  the  second  is 


horizontal  gradient  production.  Similarly,  the  steady-state  salinity  variance  should  scale  as 


-2  fh  u's'—dz  -  2  f h  w's'—dz 

S'2>  = - ^ ^ =  -  2  -<£“*-  [ h  u's'dz  (6.B.2) 

Wj<£>1/2  m1<E>1/2  J° 

Since  there  is  no  surface  salinity  flux,  the  only  contribution  using  one-dimensional 


assumptions  is  horizontal  gradient  production.  LES  turbulent  fluxes  u  'S '  and  u  /0/  were 


used  with  these  equations  to  calculate  the  values  in  Table  9. 


Temperature 

Vertical 

Horizontal 

Salinity 

Horizontal 

Variance 

Gradient 

Gradient 

Variance 

Gradient 

e'2 

6'2 

0'2 

S'2 

S'2 

Production 

Production 

Production 

gr5 

4.87e-5 

3.01e-6 

0 

1.10e-10 

0 

gr5200 

1.99e-5 

1.23e-6 

0 

1.76e-10 

0 

gr4 

5.01e-5 

2.93e-6 

-1.75e-8 

3.83e-9 

3.68e-10 

gr9 

2.08e-5 

1.19e-6 

6.29e-9 

4.50e-9 

3.59e-10 

gr6 

1.29e-3 

2.34e-6 

5.39e-5 

3.16e-6 

1.32e-7 

Table  9.  Variance  budget  contributions. 


If  the  relationship  between  temperature  variance  and  its  surface  production  in  the  no 

m3 

gradient  cases  is  used  to  estimate  —  ,  a  value  of  16.2  is  found.  This  value  gives  good 

mi 

agreement  between  surface  production  and  total  variance  for  the  weak  gradient  cases,  but  not 
for  the  strong  gradient,  which  is  not  in  steady  state.  Since  horizontal  gradients  are  the  only 
source  of  salinity  variance,  it  is  reasonable  to  estimate  a  constant  of  proportionality  between 
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total  salinity  variance  and  its  horizontal  gradient  production  as  well;  the  value  is  10.4  for  the 
400  W/m2  case  and  12.5  for  the  200  W/m2  case.  If  a  similar  constant  relates  temperature 
variance  to  its  gradient  production,  then  the  weak  gradient  case  has  variance  production 
dominated  by  surface  cooling,  and  advection  accounts  for  a  few  percent  of  the  variance. 

Fluxes  for  the  strong  gradient  case  are  shown  in  Figure  47.  Bulk  fluxes  in  time  show 
that  the  system  is  not  in  steady  state.  The  salinity  flux  profiles  have  the  same  shapes  as  the 
salinity  flux  profiles  in  the  weak  gradient  case,  but  are  more  than  an  order  of  magnitude 
larger.  Eastward  flow  through  the  upper  300m  brings  in  more  saline  water,  and  westward 
flow  below  brings  in  fresher  water,  both  giving  positive  flux  values.  Mean  geostrophic 
current  is  southward  everywhere,  with  no>>  gradient,  but  the  continual  salinization  over  most 
of  the  column,  with  freshening  near  the  bottom,  results  in  net  transport  of  salt  south.  The 
shapes  of  the  temperature  flux  profiles  resemble  the  salinity  flux  profiles  rather  than  the 
weak  gradient  case’s  temperature  profiles.  This  is  because  the  temperature  gradient  is  strong 
enough  to  dominate  over  cooling  in  determining  the  profiles  of  horizontal  flux. 

4.  Stability 

Figure  48  shows  profiles  of  detrended  temperature  and  buoyancy.  The  shapes  of  the 
temperature  and  buoyancy  profiles  remain  relatively  unchanged  by  the  gradient.  The  strong 
gradient  case  shows  that  profiles  of  temperature  and  buoyancy  are  not  in  steady  state,  as 
overturning  continues  to  change  the  shape  of  the  profiles.  Figure  49  shows  profiles  of  a  more 
sensitive  parameter,  the  buoyancy  frequency  N  squared.  All  the  simulations,  with  or  without 
horizontal  gradients,  show  an  unstable  profile  over  the  upper  250m  due  to  surface  cooling 
and  a  stable  region  near  the  bottom.  Steady-state  N 2  profiles  in  the  simulations  with 
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gradients  should  be  interpreted  to  mean  that  the  tendency  for  turbulence  to  mix  to  the  typical 
shape  without  gradients  and  the  tendency  for  overturning  to  stratify  the  system  are  in  balance. 
The  weak  gradient  cases  have  reached  this  balance,  but  in  the  strong  gradient  case  the 
overturning  is  stronger,  and  the  stratification  continues  to  increase  in  time.  The  cases  with 
400  W/m2  cooling  are  more  stable  near  the  bottom  than  are  the  cases  with  200  W/m2 
cooling.  Stronger  cooling  means  parcels  become  colder  while  at  the  surface  before  sinking. 
On  average,  there  is  a  greater  difference  between  mid-layer  parcels  and  parcels  that  have 
sunk  to  the  bottom  with  strong  cooling  than  parcels  more  weakly  cooled  at  the  surface, 
resulting  in  greater  average  stratification  near  the  bottom. 


C.  POST-CONVECTION  SPIN-DOWN 


When  cooling  at  the  surface  is  stopped,  TKE  drops  rapidly  to  an  almost  constant 
value  (Figure  50).  Turbulent  vertical  energy  becomes  less  than  ten  percent  of  the  total 
energy,  and  vertical  mixing  no  longer  opposes  die  tendency  for  overturning.  Figure  5 1  shows 
the  potential  buoyancy  field  and  horizontal  flow  at  40m  depth  48  hours  after  cooling  has 
stopped  in  the  experiment  with  the  weak  gradient.  Figure  52  is  the  same  except  it  is  for  the 
strong  gradient.  The  scale  of  the  rotating  features  is  larger  in  the  stronger  gradient,  and  may 
be  related  to  baroclinic  instability.  M2  is  defined  analogously  to  N2  as 


dx 


and  the  slope  of  a  buoyancy  surface  is 


(6.B.3) 


(6.B.4) 
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Figure  53  depicts  values  for  M2  and  N2  at  the  end  of  convection.  However,  it  is  unclear 
whether  sb  can  be  a  useful  parameter  to  describe  the  tendency  for  baroclinic  instability  when 
convection  stops  in  these  cases;  the  instability  analyses  that  use  the  slope  of  a  buoyancy 
surface  to  predict  growth  of  perturbations  assume  N2  positive  (Haine  and  Marshall,  1998). 
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VII.  CONCLUSIONS  AND  RECOMMENDATIONS 

A.  CONCLUSIONS 

One-dimensional  modeling  of  the  TKE  budget  for  the  Greenland  Sea  entraining 
mixed  layer  studied  here  agrees  well  with  large  eddy  simulation  when  universal  model 
constants  for  dissipation,  shear  production,  and  pressure  redistribution  are  tuned  for  universal 
geographically-independent  application.  The  generalized  TKE  bulk  closure  gives  a  similar 
distribution  of  TKE  components  to  those  found  in  LES  as  well  as  a  similar  ratio  of 
entrainment  buoyancy  flux  to  surface  buoyancy  flux  for  the  period  modeled.  The  ratios  are 
found  to  be  0.42  for  the  one-dimensional  model,  and  0.36  for  the  LES;  they  are  unexpectedly 
large  because  of  strong  shear  production  and  the  impact  of  thermobaricity. 

Convective  plumes  are  able  to  be  detrained  from  the  mixed  layer  even  with  stratified 
fluid  below.  A  critical  depth  hcr  and  critical  velocity  wcr  hypothesized  by  Garwood  et  al. 
(1994)  are  useful  indicators  of  the  timing  of  the  onset  of  detrainment.  The  skewness  of 
vertical  velocity  in  a  horizontal  slice  just  below  the  mixed  layer  is  a  better  parameter  than  is 
buoyancy  flux  below  the  mixed  layer  for  detection  of  detrainment  events.  The  assumptions 
about  plume  temperature  and  salinity  distributions  and  plume  vertical  velocity  used  in  the 
parameterization  of  Paluskiewicz  and  Romea  (1997)  are  justified  based  on  the  LES 
simulations.  Effects  of  detrainment  between  LES  and  this  parameterization  were  not 
compared.  The  effect  of  detrainment  on  the  fluid  between  the  mixed-layer  interface  and  the 
depth  of  neutral  buoyancy  for  the  detrained  parcel  differs  greatly  between  cases  with  zero  and 
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nonzero  stratification  below  the  mixed  layer.  Using  the  fraction  of  mixed-layer  water  present 
to  indicate  the  amount  of  plume  entrainment,  plumes  in  the  stratified  cases  entrain  much 
more. 

Convection  with  large-scale  horizontal  gradients,  which  brings  together  parcels 
originally  from  disparate  environments,  produces  greater  temperature  and  salinity  variance 
than  does  convection  with  no  large-scale  gradients.  Convection  with  horizontal  gradients 
also  increases  TKE,  mainly  in  the  horizontal  components.  The  additional  horizontal  TKE 
is  present  at  scales  larger  than  the  convective  plume  scale,  and  the  scale  for  the  horizontal 
energy  increases  with  the  strength  of  the  gradient.  Thus  the  turbulence  is  less  isotropic,  but 
the  increased  anisotropy  is  found  in  the  larger  scales.  Since  the  amount  of  vertical  TKE  is 
relatively  unaffected  by  the  presence  of  the  horizontal  gradient,  one-dimensional  TKE  and 
entrainment  models  can  safely  neglect  large-scale  horizontal  gradients  when  the  modeled 
region  is  understood  to  be  moving  with  the  flow,  and  the  water  column  is  not  subject  to 
significant  shear. 

Mean  horizontal  velocities  show  that  large-scale  overturning  occurs  simultaneously 
with  convection.  The  balance  between  the  two  processes  results  in  greater  stratification 
within  the  mixed  layer  than  in  a  mixed  layer  with  no  horizontal  gradients,  at  least  in  cases 
with  relatively  large  vertical  buoyancy  gradients.  The  trajectory  of  a  parcel  embedded  in  a 
convective  plume  during  large-scale  overturning  could  be  a  slanted  trajectory  similar  to  that 
envisioned  in  recent  literature  that  discusses  slantwise  convection.  Analysis  of  the  sources 
of  scalar  variance  in  regions  with  weak  horizontal  gradients  and  no  entrainment  shows  that 
the  variance  production  is  dominated  by  surface  processes. 
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In  the  spin-down  of  TKE  when  surface  forcing  is  stopped,  TKE  falls  to  a  steady 


value;  horizontal  gradients  affect  the  amount  of  TKE  that  remains. 


B.  RECOMMENDATIONS 


1.  One-Dimensional  Modeling 

The  one-dimensional  model  of  the  Greenland  Sea  entraining  mixed  layer  was  able 
to  produce  realistic  balance  among  TKE  components  and  among  the  TKE  budget  terms. 
However,  inclusion  of  the  tendency  for  wind  stirring  to  more  effectively  deepen  shallow 
mixed  layers  than  deeper  ones  would  be  an  important  advance  in  one-dimensional  mixed- 
layer  modeling.  Also,  the  inclusion  of  prescribed  horizontal  gradients  in  the  one-dimensional 
model  would  give  improved  variance  predictions. 

2.  Large-Scale  LES 

The  problem  of  large-scale  and  convective-scale  interaction  will  eventually  require 
a  departure  from  the  doubly-periodic  boundary  conditions  of  the  LES  model  used  in  this 
work,  especially  if  spatially-varying  rates  of  entrainment  are  to  be  allowed. 

Further  simulations  of  convection  with  large-scale  gradients  that  include  drifters 
would  be  useful  to  test  how  well  LES  fields  fit  the  descriptions  of  slantwise  convection 
found  in  recent  literature. 

Future  work  should  also  focus  on  the  horizontal  scales  produced,  what  conditions 
allow  baroclinic  instabilities  to  exist  simultaneously  with  convection,  and  how  potential 
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energy  is  converted  into  horizontal  turbulence  when  baroclinic  instability  processes  are  not 
evident. 

3.  Restratification 

Study  of  the  spin-down  of  TKE  after  surface  cooling  stops  and  its  dependence  on 
horizontal  gradients  could  improve  understanding  of  the  restratification  process. 

Simulations  of  conditions  in  which  thermobaric  interleaving  is  expected  to  occur 
would  test  the  hypotheses  that  warm  intrusions  are  preferentially  preserved,  and  have 
potential  to  explain  historical  and  current  observations. 
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APPENDIX.  FIGURES 


Figure  2.  Thermobaric  Plume  Schematic  shows  parameters  used  in  theoretical  dis¬ 
cussion.  A  cold  (near  freezing)  and  fresh  mixed  layer  overlies  a  relatively  warm  and 
saline  layer.  If  a  mixed-layer  parcel  has  downward  velocity  equal  to  or  greater  than 
the  critical  velocity  wcr,  its  momentum  will  carry  it  to  the  critical  depth  where  its 
relative  buoyancy  will  be  negative,  and  it  will  be  accelerated  downward. 
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Figure  3.  Thermobaric  Interleaving  Schematic  shows  parameters  used  in  theoretical 
discussion.  A  density-compensated  warm  and  saline  intrusion  will  be  stable  to  ver¬ 
tical  displacements.  A  cold  and  fresh  intrusion  will  be  stable  only  if  the  background 
stratification  is  stronger  than  a  critical  value  N^. 
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Figure  5.  Mixed-layer  deepening, t« 
dimensional  model  and  by  LES. 
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Figure  6.  Turbulent  kinetic  energy  (TKE)  components  as  predicted  by  the  one¬ 
dimensional  model  and  by  LES. 
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Figure  11.  TKE  contributions  as  computed  by  the  one-dimensional  model  and  by 
LES.  Differs  from  Figure  9  in  that  surface  and  entrainment  contributions  to  buoyant 
production/damping  have  been  summed. 
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Figure  12.  Initial  profiles  of  temperature  and  salinity  for  the  thermobaric  plume 
experiments. 


Relative  Buoyancy  of  a  Mixed-Layer  Parcel 


Figure  13.  Profiles  of  buoyancy  for  a  mixed-layer  parcel  if  plunged  to  any  depth  below 
the  initial  interface. 
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Figure  15.  Energy-density  spectra  for  vertical  velocity  at  95m  ten  hours  after  the 
onset  of  detrainment  for  two  unstratified  experiments.  Top:  gsea3  Bottom:  gsea3.2 
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wr  Spectrum  at  95m  -  gstrat2 


V 

Figure  17.  A  field  of  detraining  thermobaric  plumes,  rendered  using  an  isosurface  of 
temperature.  Width  of  domain  is  2.56km.  Plumes  are  descending  about  400m  below 
the  interface,  which  is  at  a  depth  of  200m. 
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Figure  18.  Time  series  from  unstratified  experiments.  Top:  LES-predicted  mixed- 
layer  depth  and  critical  depth  for  thermobaric  plume  formation.  Center:  Bulk  mixed- 
layer  Wnns  and  critical  velocity.  Bottom:  Bulk  buoyant  production  below  the  mixed 
layer. 
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Buoyant  Prod  Before  Detrainment 


Buoyant  Prod  During  Detrainment 


Figure  19.  Profiles  from  unstratified  experiments.  Bouyant  production  (buoyancy 
flux)  before  detrainment  (left)  and  during  detrainment  (right). 
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Figure  20.  Time  series  from  experiments  varying  the  stratification.  Dashed  lines 
are  critical  values.  Top:  LES-predicted  mixed-layer  depth  and  critical  depth  for 
thermobaric  plume  formation.  Center:  Bulk  mixed-layer  wrins  and  critical  velocity. 
Bottom:  Bulk  buoyant  production  below  the  mixed  layer. 
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Figure  21.  Profiles  from  experiments  varying  the  stratification.  Bouyant  production 
(buoyancy  flux)  before  detrainment  (left)  and  during  detrainment  (right). 


Skewness  of  Vertical  Velocity  at  285m 
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Figure  22.  Skewness  of  Vertical  Velocity  at  285m.  Top:  from  unstratified  experi¬ 
ments.  Bottom:  from  experiments  varying  the  stratification. 
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Figure  23.  Plots  of  LES  horizontal  mean  temperature  as  a  function  of  depth  and 
time,  from  unstratified  experiments.  Top:  Rotation  and  Full  Thermobaricity,  Center: 
No  Rotation,  Bottom:  Thermobaricity  Off  Below  the  Mixed  Layer 
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Figure  24.  Plots  of  LES  horizontal  mean  temperature  as  a  function  of  depth  and  time, 
from  experiments  varying  the  stratification.  Top:  Unstratified,  Upper  Center:  Weak 
Stratification,  Lower  Center:  Moderate  Stratification,  Bottom:  Strong  Stratification. 
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Figure  25.  Profiles  of  temperature  at  the  initial  time  and  10  hours  after  the  onset  of 
detrainment. 
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Figure  26.  Profiles  of  the  fraction  of  mixed-layer  water  present  10  and  20  hours  after 
the  onset  of  detrainment.  Solid  lines  are  at  10  hours;  dashed  lines  are  at  20  hours. 


Fraction  of  Mixed-Layer  Water  10  and  20  Hours  After  Onset  of  Detrainment 
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Figure  27.  Horizontal  cross-sections  of  temperature  (left  panel)  and  salinity  (right 
panel)  from  gsea3.2  at  mid-depth  (400m).  Box  is  2.56km  on  a  side.  Contour  inter¬ 
vals  are  .1°C  and  .006  psu  respectively.  Blue  areas  are  the  relatively  cool  and  fresh 
detraining  plumes. 


Figure  28.  Horizontal  cross-sections  of  temperature  (left  panel)  and  salinity  (left 
panel)  from  gstrat2  at  mid-depth  (400m).  Box  is  2.56km  on  a  side.  Contour  inter¬ 
vals  are  .1°C  and  .006  psu  respectively.  Blue  areas  are  the  relatively  cool  and  fresh 
detraining  plumes. 


Figure  30.  Vertical  cross-section  of  vertical  velocity  from  gstrat2.  Box  is  2.56km  wide 
by  800m  deep.  Contour  interval  is  .01  m/s.  Blue  areas  are  the  downward  velocities 
of  detraining  plumes. 
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Figure  31.  Conditions  used  in  the  horizontal  gradient  experiments.  Top:  Buoyancy 
Gradient  as  a  function  of  depth.  Bottom:  Geostrophic  Velocity  resulting  from  thermal 
wind  shear. 
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Figure  32.  Energy-density  spectra  for  vertical  velocity  at  190m  67  hours  after  cool¬ 
ing  begins.  Top:  No  Gradient  (gr5)  Center:  Weak  Gradient  (gr4)  Bottom:  Strong 
Gradient  (gr6) 


Figure  34.  LES  Salinity  at  a  depth  of  40m  from  the  weak  gradient  experiment, 
is  2.56km  on  a  side.  The  color  scale  is  (S-34.85)e-3  psu. 
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Figure  35.  LES  Vertical  Velocity  at  a  depth  of  40m  from  the  weak  gradient  experi¬ 
ment.  Box  is  2.56km  on  a  side.  The  color  scale  is  in  m/s. 
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Figure  36.  LES  Temperature  at  a  depth  of  40m  from  the  weak  gradient  experiment 
Box  is  2.56km  on  a  side.  The  color  scale  is  in  °Ce-3. 
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Bulk  Mean  Vertical  TKE 


Figure  37.  Bulk  Mean  Turbulent  Energy  as  a  function  of  time.  Top:  Total  TKE. 
Bottom:  Vertical  component  of  TKE. 
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Vertical  TKE/Total  TKE  (200  W/rrf  Cooling) 


Vertical  TKE/Total  TKE  (400  W/rrf  Cooling) 


Figure  38.  Ratio  of  Vertical  to  Total  TKE. 
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Figure  39.  Bulk  Temperature  and  Salinity  Variances  from  horizontal  means. 
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Figure  40.  LES  Salinity  at  20m  depth  from  the  strong  gradient  experiment  when  the 
box  mode  is  large.  Box  is  2.56km  on  a  side.  The  color  scale  is  in  (S-34.85)e-3  psu. 
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Figure  41.  LES  nonhydrostatic  pressure  per  unit  density  and  horizontal  flow  fields 
at  20m  depth  from  the  strong  gradient  experiment  when  the  box  mode  is  large.  Box 
is  2.56km  on  a  side.  The  color  scale  is  in  m2  s2  e-3.  The  longest  arrows  represent 
horizontal  speeds  of  9  cm/s. 
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Figure  42.  LES  Vertical  Ve' 
when  the  box  mode  is  large. 
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Figure  43.  Mean  Velocity  Profiles. 
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Figure  44.  Schematic  showing  enhancement  of  geostrophic  currents. 
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Figure  45.  Horizontal-mean  Velocity  at  the  first  grid  level. 
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Figure  46.  Horizontal  Flux  Profiles. 
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Figure  47.  Horizontal  Flux  plots  for  the  Strong  Gradient  case. 


Detrended  Temperature  Profiles  Detrended  Buoyancy  Profiles 


Figure  49.  Stability  (N2)  Profiles.  Top  is  from  the  experiments  with  200  W/m2. 
Center  is  from  the  experiments  with  400  W/m2.  Bottom  is  from  the  experiment  with 
the  strong  gradient  and  400  W/m2. 
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Figure  50.  Postconvection  Plots  of  TKE  behavior  in  time.  Top:  Bulk  Mean  TKE 
Bottom:  Vertical  TKE/Total  TKE 


lE 


a AUii///,, 


\nWv\v  I, 


BUOY 

!-508. 
-511.35 
-514.20 
-517.05 
-519.90 


Figure  51.  Potential  buoyancy  and  horizontal  flow  40m  below  the  surface  in  the  weak 
gradient  experiment  48  hours  after  cooling  stops.  Box  is  2.56km  on  a  side.  Color  scale 
is  in  m/s2  e-3.  The  longest  arrows  represent  horizontal  speeds  of  3  cm/s. 
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Figure  52.  Potential  buoyancy  and  horizontal  flow  40m  below  the  surface  in  the 
strong  gradient  experiment  48  hours  after  cooling  stops.  Box  is  2.56km  on  a  side. 
Color  scale  is  in  m/s2  e-3.  The  longest  arrows  represent  horizontal  speeds  of  3.5  cm/s. 
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Figure  53.  Postconvection  Plots  of  N2  and  M2  profiles  55  hours  after  cooling  stops 
Dashed  lines  are  bulk  averages. 
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